R Notebook
#Functions and packages
library(relaimpo)## Loading required package: MASS## Loading required package: boot## Loading required package: survey## Loading required package: grid## Loading required package: Matrix## Loading required package: survival##
## Attaching package: 'survival'## The following object is masked from 'package:boot':
##
## aml##
## Attaching package: 'survey'## The following object is masked from 'package:graphics':
##
## dotchart## Loading required package: mitools## This is the global version of package relaimpo.## If you are a non-US user, a version with the interesting additional metric pmvd is available## from Ulrike Groempings web site at prof.beuth-hochschule.de/groemping.library(ggplot2)
library(ggtext)
library(plyr)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## The following object is masked from 'package:MASS':
##
## select## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(reshape2)
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinelibrary(cowplot)
library(lme4)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:boot':
##
## logitlibrary(doBy)##
## Attaching package: 'doBy'## The following object is masked from 'package:dplyr':
##
## order_bylibrary(bipartite)## Loading required package: vegan## Loading required package: permute## Loading required package: lattice##
## Attaching package: 'lattice'## The following object is masked from 'package:boot':
##
## melanoma## This is vegan 2.5-7##
## Attaching package: 'vegan'## The following object is masked from 'package:survey':
##
## calibrate## Loading required package: sna## Loading required package: statnet.common##
## Attaching package: 'statnet.common'## The following objects are masked from 'package:base':
##
## attr, order## Loading required package: network##
## 'network' 1.17.1 (2021-06-12), part of the Statnet Project
## * 'news(package="network")' for changes since last version
## * 'citation("network")' for citation information
## * 'https://statnet.org' for help, support, and other information##
## Attaching package: 'network'## The following object is masked from 'package:plyr':
##
## is.discrete## sna: Tools for Social Network Analysis
## Version 2.6 created on 2020-10-5.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
## For citation information, type citation("sna").
## Type help(package="sna") to get started.## This is bipartite 2.16.
## For latest changes see versionlog in ?"bipartite-package". For citation see: citation("bipartite").
## Have a nice time plotting and analysing two-mode networks.##
## Attaching package: 'bipartite'## The following object is masked from 'package:vegan':
##
## nullmodel## The following object is masked from 'package:plyr':
##
## emptylibrary(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ tibble 3.1.4 ✓ purrr 0.3.4
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 2.0.1 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::arrange() masks plyr::arrange()
## x gridExtra::combine() masks dplyr::combine()
## x purrr::compact() masks plyr::compact()
## x dplyr::count() masks plyr::count()
## x tidyr::expand() masks Matrix::expand()
## x dplyr::failwith() masks plyr::failwith()
## x dplyr::filter() masks stats::filter()
## x dplyr::id() masks plyr::id()
## x dplyr::lag() masks stats::lag()
## x dplyr::mutate() masks plyr::mutate()
## x doBy::order_by() masks dplyr::order_by()
## x tidyr::pack() masks Matrix::pack()
## x car::recode() masks dplyr::recode()
## x dplyr::rename() masks plyr::rename()
## x dplyr::select() masks MASS::select()
## x purrr::some() masks car::some()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()
## x tidyr::unpack() masks Matrix::unpack()library(ggforce)
library(ggpubr)##
## Attaching package: 'ggpubr'## The following object is masked from 'package:cowplot':
##
## get_legend## The following object is masked from 'package:plyr':
##
## mutatelibrary(ggeffects)##
## Attaching package: 'ggeffects'## The following object is masked from 'package:cowplot':
##
## get_titlelibrary(viridis)## Loading required package: viridisLitelibrary(viridisLite)
library(phia)
library(jtools)
library(EcoSimR)
library(cooccur)
library(fitdistrplus)
theme_ines<-theme(axis.text = element_text(size=14), axis.title = element_text(size=14, face="bold"), legend.text = element_text(size=12), strip.text = element_text(size=14), plot.title = element_text(size=14, face="bold"), panel.grid=element_line(colour="white"), panel.background = element_rect(fill="white") , axis.line = element_line(size = 0.5, linetype = "solid",
colour = "black"), strip.background = element_rect(fill="white"))
save_plot<-function(dir, width=15, height=10, ...){
ggsave(dir, width = width, height = height, units = c("cm"))
}1 - Importing data
Importing coexistence data and the leaf age
df<-read.csv("./Data/Coexistence_data_final_corr.csv", header=TRUE, row.names = NULL, as.is=c(2))
dir.create("./Plots")## Warning in dir.create("./Plots"): './Plots' already existsleaf_conversion<-read.table("./Data/leaf_conversion.txt", header=TRUE, row.names = NULL, as.is=c(1,4))
df$leaf_age<-sapply(c(1:length(df[,1])), function(x){
leaf_conversion$Age[which(as.character(leaf_conversion$Treatment)==as.character(df$Treatment[x]) & leaf_conversion$Box==df$Box[x] & leaf_conversion$Leaf==df$Leaf[x])]
})
df$Density<-mapvalues(df$Treatment, as.character(levels(as.factor(df$Treatment))), c("10:10","10:10", "10:10", "19:1", "19:1", "1:19","1:19", "20:0", "0:20"))
df$Timing<-mapvalues(df$Treatment, as.character(levels(as.factor(df$Treatment))), c("same time", "Tu first", "Te first", "same time", "Te first", "same time", "Tu first", "single", "single"))
str(df)## 'data.frame': 3214 obs. of 11 variables:
## $ Box : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment: chr "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
## $ Block : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Leaf : int 2 2 2 2 2 2 3 3 3 3 ...
## $ Leaflet : int 1 2 3 1 2 3 1 2 3 4 ...
## $ Direction: Factor w/ 2 levels "Down","Up": 1 1 1 2 2 2 1 1 1 1 ...
## $ Te : int 4 2 1 6 6 3 0 0 4 0 ...
## $ Tu : int 1 2 5 1 12 6 1 2 0 2 ...
## $ leaf_age : chr "old" "old" "old" "old" ...
## $ Density : chr "10:10" "10:10" "10:10" "10:10" ...
## $ Timing : chr "same time" "same time" "same time" "same time" ...Calculating proportions and ratios
df$ratio<-sapply(c(1:length(df[,1])), function(x){
#print(x)
if(df$Treatment[x]=="20Te" | df$Treatment[x]=="20Tu"){
a<-NA
}else{
if(df$Te[x]==0){
te<-1
}else{
te<-df$Te[x]
}
if(df$Tu[x]==0){
tu<-1
}else{
tu<-df$Tu[x]
}
a<-te/tu
}
a
})
df$Total<-sapply(c(1:length(df[,1])), function(x){
if(df$Treatment[x]=="20Te"){
a<-df$Te[x]
}else if(df$Treatment[x]=="20Tu"){
a<-df$Tu[x]
}else{
a<-df$Te[x]+df$Tu[x]
}
a
})
df$proportionTe<-sapply(c(1:length(df[,1])), function(x){
if(df$Treatment[x]=="20Tu"){
a<-0
}else if(df$Te[x]==0 |df$Total[x]==0 ){
a<-0
} else{
te<-df$Te[x]
a<-te/df$Total[x]
}
a
})
df$proportionTu<-sapply(c(1:length(df[,1])), function(x){
if(df$Treatment[x]=="20Te"){
a<-0
}else if(df$Tu[x]==0 |df$Total[x]==0 ){
a<-0
} else{
tu<-df$Tu[x]
a<-tu/df$Total[x]
}
a
})Setting up data frame with data per box
For the statistical analysis we will analyse data per box first. Here we are summing the number of individuals of each species per box - df_box To test a different distribution per leaf we will also create a data frame with number of individuals per leaf and per direction - df_box_leaf_dir
df_box_leaf_dir<-df %>%
group_by(Treatment, Density, Timing , Block, Box, Direction,Leaf, leaf_age) %>%
summarise(Te=sum(Te, na.rm=TRUE), Tu=sum(Tu, na.rm=TRUE)) %>%
as.data.frame()## `summarise()` has grouped output by 'Treatment', 'Density', 'Timing', 'Block', 'Box', 'Direction', 'Leaf'. You can override using the `.groups` argument.df_box<-df %>%
group_by(Treatment,Density, Timing , Block, Box) %>%
summarise(Te=sum(Te, na.rm=TRUE), Tu=sum(Tu, na.rm=TRUE)) %>%
as.data.frame()## `summarise()` has grouped output by 'Treatment', 'Density', 'Timing', 'Block'. You can override using the `.groups` argument.head(df_box_leaf_dir)## Lets add our usual proportion information
df_box$ratio<-sapply(c(1:length(df_box[,1])), function(x){
#print(x)
if(df_box$Treatment[x]=="20Te" | df_box$Treatment[x]=="20Tu"){
a<-NA
}else{
if(df_box$Te[x]==0){
te<-1
}else{
te<-df_box$Te[x]
}
if(df_box$Tu[x]==0){
tu<-1
}else{
tu<-df_box$Tu[x]
}
a<-te/tu
}
a
})
df_box_leaf_dir$ratio<-sapply(c(1:length(df_box_leaf_dir[,1])), function(x){
#print(x)
if(df_box_leaf_dir$Treatment[x]=="20Te" | df_box_leaf_dir$Treatment[x]=="20Tu"){
a<-NA
}else{
if(df_box_leaf_dir$Te[x]==0){
te<-1
}else{
te<-df_box_leaf_dir$Te[x]
}
if(df_box_leaf_dir$Tu[x]==0){
tu<-1
}else{
tu<-df_box_leaf_dir$Tu[x]
}
a<-te/tu
}
a
})
df_box_leaf_dir$Total<-sapply(c(1:length(df_box_leaf_dir[,1])), function(x){
if(df_box_leaf_dir$Treatment[x]=="20Te"){
a<-df_box_leaf_dir$Te[x]
}else if(df_box_leaf_dir$Treatment[x]=="20Tu"){
a<-df_box_leaf_dir$Tu[x]
}else{
a<-df_box_leaf_dir$Te[x]+df_box_leaf_dir$Tu[x]
}
a
})
df_box_leaf_dir$proportionTe<-sapply(c(1:length(df_box_leaf_dir[,1])), function(x){
if(df_box_leaf_dir$Treatment[x]=="20Tu"){
a<-0
}else if(df_box_leaf_dir$Te[x]==0 |df_box_leaf_dir$Total[x]==0 ){
a<-0
} else{
te<-df_box_leaf_dir$Te[x]
a<-te/df_box_leaf_dir$Total[x]
}
a
})
df_box_leaf_dir$proportionTu<-sapply(c(1:length(df_box_leaf_dir[,1])), function(x){
if(df_box_leaf_dir$Treatment[x]=="20Te"){
a<-0
}else if(df_box_leaf_dir$Tu[x]==0 |df_box_leaf_dir$Total[x]==0 ){
a<-0
} else{
tu<-df_box_leaf_dir$Tu[x]
a<-tu/df_box_leaf_dir$Total[x]
}
a
})
df_box$Total<-sapply(c(1:length(df_box[,1])), function(x){
if(df_box$Treatment[x]=="20Te"){
a<-df_box$Te[x]
}else if(df_box$Treatment[x]=="20Tu"){
a<-df_box$Tu[x]
}else{
a<-df_box$Te[x]+df_box$Tu[x]
}
a
})
df_box$proportionTe<-sapply(c(1:length(df_box[,1])), function(x){
if(df_box$Treatment[x]=="20Tu"){
a<-0
}else if(df_box$Te[x]==0 |df_box$Total[x]==0 ){
a<-0
} else{
te<-df_box$Te[x]
a<-te/df_box$Total[x]
}
a
})
df_box$proportionTu<-sapply(c(1:length(df_box[,1])), function(x){
if(df_box$Treatment[x]=="20Te"){
a<-0
}else if(df_box$Tu[x]==0 |df_box$Total[x]==0 ){
a<-0
} else{
tu<-df_box$Tu[x]
a<-tu/df_box$Total[x]
}
a
})Checking distribution of data
The variables that we are interested to know the distribution are: Te and Tu (because of the binomial analyses)
hist(df$Te)
hist(df$Tu)
descdist(df$Te[-which(is.na(df$Te/df$Total))], discrete = FALSE, boot=500)
## summary statistics
## ------
## min: 0 max: 149
## median: 4
## mean: 10.08053
## estimated sd: 15.99749
## estimated skewness: 3.468693
## estimated kurtosis: 19.60968descdist(df$Tu[-which(is.na(df$Tu/df$Total))], discrete = FALSE, boot=500)
## summary statistics
## ------
## min: 0 max: 125
## median: 2
## mean: 4.420565
## estimated sd: 8.373305
## estimated skewness: 5.241611
## estimated kurtosis: 47.52895Clearly not a normal distribution so we will use the binomial family in glmer to analyse data
Statistical analyses
2 - Does the timing and initial density change the probability of coexistence?
2.1 - Whole data set
Figure 1
df_box$Density2<-factor(df_box$Density, levels=c("1:19","10:10","19:1"))
ggplot(subset(df_box, Treatment!="20Te" & Treatment!="20Tu"), aes(x=Density2, y=proportionTe, fill=interaction(Density2,Timing)))+
facet_wrap("Timing", ncol=3, scales = "free_x")+
geom_hline(yintercept = 0.5)+
geom_boxplot()+
theme_ines+
scale_fill_manual(values = c("#e0f3f8", "#91bfdb", "#4575b4","#d73027", "#fc8d59", "#fee090", "#ffffbf"), name="Treatment", guide="none", drop=TRUE)+
ylab(c("Proportion of T. evansi"))+
xlab(c("Initial T. evansi: T. urticae density"))+
#scale_x_discrete(labels=c())+
theme(legend.title = element_text(size=18, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75))+
stat_summary(fun.y=median, geom="smooth", aes(group=Timing, colour=Timing), lwd=1)+
scale_colour_manual(values = c( "#4575b4","#d73027", "#fee090"), name="Treatment", guide="none")## Warning: `fun.y` is deprecated. Use `fun` instead.
save_plot("./Plots/Figure1.pdf", width=15, height=10)Stats
Random part of the model
First we need to remove the single species regimes, to test for the effects Since here we only have as random variable Block, our model will include block as random factor Our unit of analysis will be the box
df$Block2<-as.factor(df$Block)
df$Box2<-as.factor(df$Box)
str(df)## 'data.frame': 3214 obs. of 17 variables:
## $ Box : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : chr "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
## $ Block : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Leaf : int 2 2 2 2 2 2 3 3 3 3 ...
## $ Leaflet : int 1 2 3 1 2 3 1 2 3 4 ...
## $ Direction : Factor w/ 2 levels "Down","Up": 1 1 1 2 2 2 1 1 1 1 ...
## $ Te : int 4 2 1 6 6 3 0 0 4 0 ...
## $ Tu : int 1 2 5 1 12 6 1 2 0 2 ...
## $ leaf_age : chr "old" "old" "old" "old" ...
## $ Density : chr "10:10" "10:10" "10:10" "10:10" ...
## $ Timing : chr "same time" "same time" "same time" "same time" ...
## $ ratio : num 4 1 0.2 6 0.5 0.5 1 0.5 4 0.5 ...
## $ Total : int 5 4 6 7 18 9 1 2 4 2 ...
## $ proportionTe: num 0.8 0.5 0.167 0.857 0.333 ...
## $ proportionTu: num 0.2 0.5 0.833 0.143 0.667 ...
## $ Block2 : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ Box2 : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...df_wsingle<-subset(df, Treatment!="20Te" & Treatment!="20Tu")
df_wsingle<-droplevels(df_wsingle)
str(df_wsingle)## 'data.frame': 2500 obs. of 17 variables:
## $ Box : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : chr "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
## $ Block : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Leaf : int 2 2 2 2 2 2 3 3 3 3 ...
## $ Leaflet : int 1 2 3 1 2 3 1 2 3 4 ...
## $ Direction : Factor w/ 2 levels "Down","Up": 1 1 1 2 2 2 1 1 1 1 ...
## $ Te : int 4 2 1 6 6 3 0 0 4 0 ...
## $ Tu : int 1 2 5 1 12 6 1 2 0 2 ...
## $ leaf_age : chr "old" "old" "old" "old" ...
## $ Density : chr "10:10" "10:10" "10:10" "10:10" ...
## $ Timing : chr "same time" "same time" "same time" "same time" ...
## $ ratio : num 4 1 0.2 6 0.5 0.5 1 0.5 4 0.5 ...
## $ Total : int 5 4 6 7 18 9 1 2 4 2 ...
## $ proportionTe: num 0.8 0.5 0.167 0.857 0.333 ...
## $ proportionTu: num 0.2 0.5 0.833 0.143 0.667 ...
## $ Block2 : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ Box2 : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...levels(as.factor(df_wsingle$Treatment))## [1] "10Te10Tu" "10Te10Tu48h" "10Te48h10Tu" "19Te1Tu" "19Te48h1Tu"
## [6] "1Te19Tu" "1Te19Tu48h"r1<-glmer(cbind(Te,Total)~1+ (1|Block2), data=df_wsingle, family = binomial(link="logit"))
summary(r1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(Te, Total) ~ 1 + (1 | Block2)
## Data: df_wsingle
##
## AIC BIC logLik deviance df.resid
## 11424.3 11436.0 -5710.2 11420.3 2498
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.0545 -0.8376 0.0000 0.5633 3.0748
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block2 (Intercept) 0.0002307 0.01519
## Number of obs: 2500, groups: Block2, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3659 0.0143 -25.58 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Complete model with Block as random factor
b1<-glmer(cbind(Te,Total)~ Density*Timing+(1|Block2), data=df_wsingle, family = binomial(link="logit"))## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients## boundary (singular) fit: see ?isSingularb2<-glmer(cbind(Te,Total)~ Density+Timing+(1|Block2), data=df_wsingle, family = binomial(link="logit"))## boundary (singular) fit: see ?isSingularb3<-glmer(cbind(Te,Total)~ Timing+(1|Block2), data=df_wsingle, family = binomial(link="logit"))
b4<-glmer(cbind(Te,Total)~ Density+(1|Block2), data=df_wsingle, family = binomial(link="logit"))## boundary (singular) fit: see ?isSingularb5<-glmer(cbind(Te,Total)~ Treatment+(1|Block2), data=df_wsingle, family = binomial(link="logit"))## boundary (singular) fit: see ?isSingularanova(b1, b2, b3, b4, b5)The lowest AIC is from models b1 and b5, which correspond to the complete model, so lets take a look at the summary of the model
summary(b1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(Te, Total) ~ Density * Timing + (1 | Block2)
## Data: df_wsingle
##
## AIC BIC logLik deviance df.resid
## 7602.9 7649.4 -3793.4 7586.9 2492
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.3247 -0.3826 0.0000 0.1353 9.5786
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block2 (Intercept) 2.142e-17 4.628e-09
## Number of obs: 2500, groups: Block2, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.30669 0.05790 -22.567 < 2e-16 ***
## Density10:10 0.81807 0.06352 12.880 < 2e-16 ***
## Density19:1 1.28768 0.06087 21.155 < 2e-16 ***
## TimingTe first -0.01084 0.02572 -0.421 0.673
## TimingTu first -1.71343 0.09913 -17.285 < 2e-16 ***
## Density10:10:TimingTe first 0.30846 0.04311 7.155 8.37e-13 ***
## Density10:10:TimingTu first 1.35317 0.10696 12.651 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) D10:10 Dn19:1 TimingTefirst TimingTufirst
## Densty10:10 -0.912
## Density19:1 -0.951 0.867
## TimingTefirst 0.000 0.000 -0.225
## TimingTufirst -0.584 0.533 0.556 0.000
## Density10:10:TimingTefirst 0.000 -0.249 0.134 -0.597 0.000
## Density10:10:TimingTufirst 0.541 -0.594 -0.515 0.000 -0.927
## Density10:10:TimingTefirst
## Densty10:10
## Density19:1
## TimingTefirst
## TimingTufirst
## Density10:10:TimingTefirst
## Density10:10:TimingTufirst 0.148
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularsummary(b5)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(Te, Total) ~ Treatment + (1 | Block2)
## Data: df_wsingle
##
## AIC BIC logLik deviance df.resid
## 7602.9 7649.4 -3793.4 7586.9 2492
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.3247 -0.3826 0.0000 0.1353 9.5786
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block2 (Intercept) 0 0
## Number of obs: 2500, groups: Block2, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.48862 0.02610 -18.719 <2e-16 ***
## Treatment10Te10Tu48h -0.36026 0.04017 -8.967 <2e-16 ***
## Treatment10Te48h10Tu 0.29763 0.03460 8.602 <2e-16 ***
## Treatment19Te1Tu 0.46961 0.03215 14.608 <2e-16 ***
## Treatment19Te48h1Tu 0.45877 0.03148 14.575 <2e-16 ***
## Treatment1Te19Tu -0.81807 0.06352 -12.880 <2e-16 ***
## Treatment1Te19Tu48h -2.53150 0.08459 -29.927 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) T10T10 T10T48 T19T1T T19T48 Tr1T19T
## Trt10T10T48 -0.650
## Trt10T4810T -0.754 0.490
## Trtmnt19T1T -0.812 0.528 0.613
## Trtm19T481T -0.829 0.539 0.626 0.673
## Trtmnt1T19T -0.411 0.267 0.310 0.334 0.341
## Trtm1T19T48 -0.309 0.200 0.233 0.251 0.256 0.127
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularAnova(b1, type=3)Anova(b5, type=3)There is a clear effect of treatment in the proportion of Te females. To test the impact of orde of arrival and initial frequency we can perform contrasts using the phia package.
Contrasts
Using a single factor
b5_2<-glmer(cbind(Te,Total)~ Treatment+(1|Block2), data=df_wsingle, family = binomial(link="logit"))## boundary (singular) fit: see ?isSingularsummary(b5_2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(Te, Total) ~ Treatment + (1 | Block2)
## Data: df_wsingle
##
## AIC BIC logLik deviance df.resid
## 7602.9 7649.4 -3793.4 7586.9 2492
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.3247 -0.3826 0.0000 0.1353 9.5786
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block2 (Intercept) 0 0
## Number of obs: 2500, groups: Block2, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.48862 0.02610 -18.719 <2e-16 ***
## Treatment10Te10Tu48h -0.36026 0.04017 -8.967 <2e-16 ***
## Treatment10Te48h10Tu 0.29763 0.03460 8.602 <2e-16 ***
## Treatment19Te1Tu 0.46961 0.03215 14.608 <2e-16 ***
## Treatment19Te48h1Tu 0.45877 0.03148 14.575 <2e-16 ***
## Treatment1Te19Tu -0.81807 0.06352 -12.880 <2e-16 ***
## Treatment1Te19Tu48h -2.53150 0.08459 -29.927 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) T10T10 T10T48 T19T1T T19T48 Tr1T19T
## Trt10T10T48 -0.650
## Trt10T4810T -0.754 0.490
## Trtmnt19T1T -0.812 0.528 0.613
## Trtm19T481T -0.829 0.539 0.626 0.673
## Trtmnt1T19T -0.411 0.267 0.310 0.334 0.341
## Trtm1T19T48 -0.309 0.200 0.233 0.251 0.256 0.127
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularlevels(as.factor(df_wsingle$Treatment))## [1] "10Te10Tu" "10Te10Tu48h" "10Te48h10Tu" "19Te1Tu" "19Te48h1Tu"
## [6] "1Te19Tu" "1Te19Tu48h"test_int2<-testInteractions(b5_2, pairwise="Treatment", adjustment = "fdr")
test_int2#Now we just select the ones we want to compare
#effect of timing
test_int2[c(1, 2,7,16, 21 ),]#effect of density
test_int2[c(3,5,11, 13,17 ),]# Testing overall order of arrival and overall density
mat_overall<-matrix(ncol=4, nrow=7)
colnames(mat_overall)<-c("Order of arrival Te first", "Order of arrival Tu first", "Te higher initially", "Tu higher initially")
rownames(mat_overall)<-c("10Te10Tu","10Te10Tu48h", "10Te48h10Tu", "19Te1Tu", "19Te48h1Tu", "1Te19Tu", "1Te19Tu48h" )
mat_overall[1,1]<- 1
mat_overall[2,1]<- 0
mat_overall[3,1]<- -1
mat_overall[4,1]<- 1
mat_overall[5,1]<- -1
mat_overall[6,1]<- 0
mat_overall[7,1]<- 0
mat_overall[1,2]<- 1
mat_overall[2,2]<- -1
mat_overall[3,2]<- 0
mat_overall[4,2]<- 0
mat_overall[5,2]<- 0
mat_overall[6,2]<- 1
mat_overall[7,2]<- -1
mat_overall[1,3]<- 1
mat_overall[2,3]<- 1
mat_overall[3,3]<- 1
mat_overall[4,3]<- -1
mat_overall[5,3]<- -1
mat_overall[6,3]<- 0
mat_overall[7,3]<- 0
mat_overall[1,4]<- 1
mat_overall[2,4]<- 1
mat_overall[3,4]<- 1
mat_overall[4,4]<- 0
mat_overall[5,4]<- 0
mat_overall[6,4]<- -1
mat_overall[7,4]<- -1
mat_overall## Order of arrival Te first Order of arrival Tu first
## 10Te10Tu 1 1
## 10Te10Tu48h 0 -1
## 10Te48h10Tu -1 0
## 19Te1Tu 1 0
## 19Te48h1Tu -1 0
## 1Te19Tu 0 1
## 1Te19Tu48h 0 -1
## Te higher initially Tu higher initially
## 10Te10Tu 1 1
## 10Te10Tu48h 1 1
## 10Te48h10Tu 1 1
## 19Te1Tu -1 0
## 19Te48h1Tu -1 0
## 1Te19Tu 0 -1
## 1Te19Tu48h 0 -1test_int3<-testInteractions(b5_2, custom=list(Treatment=mat_overall), adjustment = "fdr")
test_int32.2 - Test the importance of order of leaves
To test this we will run the model for the subset of the data with leaves 2 -4 first or 3 - 5 first
str(df_wsingle)## 'data.frame': 2500 obs. of 17 variables:
## $ Box : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : chr "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
## $ Block : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Leaf : int 2 2 2 2 2 2 3 3 3 3 ...
## $ Leaflet : int 1 2 3 1 2 3 1 2 3 4 ...
## $ Direction : Factor w/ 2 levels "Down","Up": 1 1 1 2 2 2 1 1 1 1 ...
## $ Te : int 4 2 1 6 6 3 0 0 4 0 ...
## $ Tu : int 1 2 5 1 12 6 1 2 0 2 ...
## $ leaf_age : chr "old" "old" "old" "old" ...
## $ Density : chr "10:10" "10:10" "10:10" "10:10" ...
## $ Timing : chr "same time" "same time" "same time" "same time" ...
## $ ratio : num 4 1 0.2 6 0.5 0.5 1 0.5 4 0.5 ...
## $ Total : int 5 4 6 7 18 9 1 2 4 2 ...
## $ proportionTe: num 0.8 0.5 0.167 0.857 0.333 ...
## $ proportionTu: num 0.2 0.5 0.833 0.143 0.667 ...
## $ Block2 : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ Box2 : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...pair2_4_old<-subset(df_wsingle, (Leaf==2 & leaf_age=="old") | (Leaf==4 & leaf_age=="old") | (Leaf==3 & leaf_age=="new") | (Leaf==5 & leaf_age=="new"))
pair2_4_new<-subset(df_wsingle, (Leaf==3 & leaf_age=="old") | (Leaf==5 & leaf_age=="old") | (Leaf==2 & leaf_age=="new") | (Leaf==4 & leaf_age=="new"))
g24Old<-glmer(cbind(Te,Total)~ Treatment+(1|Block2), data=pair2_4_old, family = binomial(link="logit"))## boundary (singular) fit: see ?isSingularsummary(g24Old)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(Te, Total) ~ Treatment + (1 | Block2)
## Data: pair2_4_old
##
## AIC BIC logLik deviance df.resid
## 3768.4 3809.5 -1876.2 3752.4 1242
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4200 -0.4321 0.0000 0.1346 9.7348
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block2 (Intercept) 0 0
## Number of obs: 1250, groups: Block2, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.52429 0.03926 -13.355 < 2e-16 ***
## Treatment10Te10Tu48h -0.11894 0.05602 -2.123 0.0337 *
## Treatment10Te48h10Tu 0.32034 0.05050 6.343 2.25e-10 ***
## Treatment19Te1Tu 0.50774 0.04666 10.881 < 2e-16 ***
## Treatment19Te48h1Tu 0.49115 0.04790 10.254 < 2e-16 ***
## Treatment1Te19Tu -0.86836 0.09743 -8.913 < 2e-16 ***
## Treatment1Te19Tu48h -2.52471 0.11774 -21.442 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) T10T10 T10T48 T19T1T T19T48 Tr1T19T
## Trt10T10T48 -0.701
## Trt10T4810T -0.777 0.545
## Trtmnt19T1T -0.841 0.590 0.654
## Trtm19T481T -0.820 0.574 0.637 0.689
## Trtmnt1T19T -0.403 0.282 0.313 0.339 0.330
## Trtm1T19T48 -0.333 0.234 0.259 0.280 0.273 0.134
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularAnova(g24Old, type=3)g24New<-glmer(cbind(Te,Total)~ Treatment+(1|Block2), data=pair2_4_new, family = binomial(link="logit"))## boundary (singular) fit: see ?isSingularsummary(g24New)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(Te, Total) ~ Treatment + (1 | Block2)
## Data: pair2_4_new
##
## AIC BIC logLik deviance df.resid
## 3788.1 3829.2 -1886.1 3772.1 1242
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.2080 -0.3468 0.0000 0.1375 5.9739
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block2 (Intercept) 8.736e-18 2.956e-09
## Number of obs: 1250, groups: Block2, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.46012 0.03496 -13.163 < 2e-16 ***
## Treatment10Te10Tu48h -0.65741 0.05938 -11.071 < 2e-16 ***
## Treatment10Te48h10Tu 0.28270 0.04772 5.924 3.14e-09 ***
## Treatment19Te1Tu 0.43806 0.04483 9.771 < 2e-16 ***
## Treatment19Te48h1Tu 0.43257 0.04180 10.350 < 2e-16 ***
## Treatment1Te19Tu -0.78099 0.08384 -9.316 < 2e-16 ***
## Treatment1Te19Tu48h -2.52713 0.12192 -20.728 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) T10T10 T10T48 T19T1T T19T48 Tr1T19T
## Trt10T10T48 -0.589
## Trt10T4810T -0.733 0.431
## Trtmnt19T1T -0.780 0.459 0.571
## Trtm19T481T -0.836 0.492 0.613 0.652
## Trtmnt1T19T -0.417 0.245 0.305 0.325 0.349
## Trtm1T19T48 -0.287 0.169 0.210 0.224 0.240 0.120
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularAnova(g24New, type=3)test_intOld<-testInteractions(g24Old, pairwise="Treatment", adjustment = "fdr")
test_intOldtest_intNew<-testInteractions(g24New, pairwise="Treatment", adjustment = "fdr")
test_intNew#Now we just select the ones we want to compare
#effect of timing
test_intOld[c(1, 2,16, 21, 7 ),]test_intNew[c(1, 2,16, 21, 7 ),]#effect of density
test_intOld[c(3,5,11, 13,17 ),]test_intNew[c(3,5,11, 13,17 ),]test_intOld2<-testInteractions(g24Old, custom=list(Treatment=mat_overall), adjustment = "fdr")
test_intNew2<-testInteractions(g24New, custom=list(Treatment=mat_overall), adjustment = "fdr")
test_intOld2test_intNew2We have the same results as before! Which indicates that adding the different pairs of leaves first does not affect the impact of order of arrival or initial frequency.
3 - Impact of spatial variation
3.1 - Is spatial variation mediating coexistence
Analysing differences to the occupancy of single treatments
This analyses do not take into account leaf age!
Here we will calculate the percentage of individuals from each box that are in each leaf for the single regime. Then we will don the same for the treatments and will then compare the difference
Fig 3, S3
# Start with calculating total number of individuals per leaf
single_df<-subset(df, Treatment=="20Te" | Treatment=="20Tu")
str(single_df)## 'data.frame': 714 obs. of 17 variables:
## $ Box : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : chr "20Te" "20Te" "20Te" "20Te" ...
## $ Block : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Leaf : int 2 2 2 2 2 2 3 3 3 3 ...
## $ Leaflet : int 1 2 3 1 2 3 1 2 3 4 ...
## $ Direction : Factor w/ 2 levels "Down","Up": 1 1 1 2 2 2 1 1 1 1 ...
## $ Te : int 2 20 0 13 10 0 7 15 2 4 ...
## $ Tu : int NA NA NA NA NA NA NA NA NA NA ...
## $ leaf_age : chr "old" "old" "old" "old" ...
## $ Density : chr "20:0" "20:0" "20:0" "20:0" ...
## $ Timing : chr "single" "single" "single" "single" ...
## $ ratio : num NA NA NA NA NA NA NA NA NA NA ...
## $ Total : int 2 20 0 13 10 0 7 15 2 4 ...
## $ proportionTe: num 1 1 0 1 1 0 1 1 1 1 ...
## $ proportionTu: num 0 0 0 0 0 0 0 0 0 0 ...
## $ Block2 : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ Box2 : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...total_box_single_leaf<-single_df %>%
group_by(Treatment, Box, Leaf, Block) %>%
summarize(sumTe=sum(Te, na.rm=TRUE),sumTu=sum(Tu, na.rm=TRUE))%>%
as.data.frame()## `summarise()` has grouped output by 'Treatment', 'Box', 'Leaf'. You can override using the `.groups` argument.# Same for the treatments
total_box_treat_leaf<-df %>%
group_by(Treatment, Box, Leaf, Block) %>%
summarize(sumTe=sum(Te, na.rm=TRUE),sumTu=sum(Tu, na.rm=TRUE))%>%
as.data.frame()## `summarise()` has grouped output by 'Treatment', 'Box', 'Leaf'. You can override using the `.groups` argument.## Calculating percentage
#For single
total_box_single_leaf$PercentTe<-sapply(c(1:length(total_box_single_leaf[,1])), function(x){
aux<-subset(total_box_single_leaf, Treatment==total_box_single_leaf$Treatment[x] & Box==total_box_single_leaf$Box[x])
aux_sum<-sum(aux$sumTe[1:4], na.rm=TRUE)
a<-total_box_single_leaf$sumTe[x]/aux_sum
if(a=="NaN"){
a<-NA
}
a
})
total_box_single_leaf$PercentTu<-sapply(c(1:length(total_box_single_leaf[,1])), function(x){
aux<-subset(total_box_single_leaf, Treatment==total_box_single_leaf$Treatment[x] & Box==total_box_single_leaf$Box[x])
aux_sum<-sum(aux$sumTu[1:4], na.rm=TRUE)
a<-total_box_single_leaf$sumTu[x]/aux_sum
if(a=="NaN"){
a<-NA
}
a
})
# For treatments
total_box_treat_leaf$PercentTe<-sapply(c(1:length(total_box_treat_leaf[,1])), function(x){
aux<-subset(total_box_treat_leaf, Treatment==total_box_treat_leaf$Treatment[x] & Box==total_box_treat_leaf$Box[x])
aux_sum<-sum(aux$sumTe[1:4], na.rm=TRUE)
a<-total_box_treat_leaf$sumTe[x]/aux_sum
if(a=="NaN"){
a<-NA
}
a
})
total_box_treat_leaf$PercentTu<-sapply(c(1:length(total_box_treat_leaf[,1])), function(x){
aux<-subset(total_box_treat_leaf, Treatment==total_box_treat_leaf$Treatment[x] & Box==total_box_treat_leaf$Box[x])
aux_sum<-sum(aux$sumTu[1:4], na.rm=TRUE)
a<-total_box_treat_leaf$sumTu[x]/aux_sum
if(a=="NaN"){
a<-NA
}
a
})
total_box_treat_leaf### Now calculating the difference to the average percentage of the single regime
# first creating average data base for the single regimes
mean_percentage_single<-total_box_single_leaf %>%
group_by(Treatment, Leaf)%>%
summarize(meanTe=mean(PercentTe), meanTu=mean(PercentTu))%>%
as.data.frame()## `summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.mean_percentage_single#Now adding to thr percent data frame
aux_diffTe<-t(sapply(c(1:length(total_box_treat_leaf[,1])), function(x){
aux<-subset(mean_percentage_single, Leaf==total_box_treat_leaf$Leaf[x])
te<-total_box_treat_leaf$PercentTe[x]-aux$meanTe[1]
tu<-total_box_treat_leaf$PercentTu[x]-aux$meanTu[2] #because the first is always NA
c(te, tu)
}))
colnames(aux_diffTe)<-c("diffTe", "diffTu")
total_box_treat_leaf<-as.data.frame(cbind(total_box_treat_leaf[,c(1:8)], aux_diffTe))
str(total_box_treat_leaf)## 'data.frame': 360 obs. of 10 variables:
## $ Treatment: chr "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
## $ Box : int 1 1 1 1 2 2 2 2 3 3 ...
## $ Leaf : int 2 3 4 5 2 3 4 5 2 3 ...
## $ Block : int 1 1 1 1 1 1 1 1 1 1 ...
## $ sumTe : int 22 22 80 8 68 39 102 28 19 22 ...
## $ sumTu : int 27 51 27 23 30 34 54 1 24 32 ...
## $ PercentTe: num 0.1667 0.1667 0.6061 0.0606 0.2869 ...
## $ PercentTu: num 0.211 0.398 0.211 0.18 0.252 ...
## $ diffTe : num 0.0459 -0.1074 0.155 -0.0935 0.1662 ...
## $ diffTu : num 0.0436 0.1188 -0.0968 -0.0656 0.0848 ...total_box_treat_leaf$Treatment2<-factor(total_box_treat_leaf$Treatment, levels=c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"))
total_box_treat_leaf$Treatment3<-mapvalues(total_box_treat_leaf$Treatment,c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"), c("1:19\nsame time", "10:10\nsame time", "19:1\nsame time", "1:19\nTu first", "10:10\nTu first","10:10\nTe first", "19:1\nTe first") )
total_box_treat_leaf$Treatment4<-factor(total_box_treat_leaf$Treatment3, levels=c("1:19\nsame time", "10:10\nsame time", "19:1\nsame time", "1:19\nTu first", "10:10\nTu first","10:10\nTe first", "19:1\nTe first") )
plot_percent_Te<-ggplot(subset(total_box_treat_leaf,Treatment!="20Te" & Treatment!="20Tu"), aes(x=as.factor(Leaf), y=diffTe, fill=Treatment4))+
facet_grid(.~Treatment4)+
geom_hline(yintercept = 0)+
geom_boxplot()+
theme_ines+
scale_fill_manual(values = c("#e0f3f8", "#91bfdb", "#4575b4", "#fee090", "#ffffbf","#d73027", "#fc8d59"), name="Treatment")+
ylab(c("Percentage of T. evansi per leaf\nrelative to single control"))+
xlab(c("Leaf Position"))+
#scale_x_discrete(labels=c())+
theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
plot_percent_Te
save_plot("./Plots/Percent_Single_Te.pdf", width=25, height=15)
plot_percent_Tu<-ggplot(subset(total_box_treat_leaf,Treatment!="20Te" & Treatment!="20Tu"), aes(x=as.factor(Leaf), y=diffTu, fill=Treatment4))+
facet_grid(.~Treatment4)+
geom_hline(yintercept = 0)+
geom_boxplot()+
theme_ines+
scale_fill_manual(values = c("#e0f3f8", "#91bfdb", "#4575b4", "#fee090", "#ffffbf","#d73027", "#fc8d59"), name="Treatment")+
ylab(c("Density of T. urticae per leaf\nrelative to single control"))+
xlab(c("Leaf Position"))+
#scale_x_discrete(labels=c())+
theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
plot_percent_Tu## Warning: Removed 4 rows containing non-finite values (stat_boxplot).
save_plot("./Plots/Percent_Single_Tu.pdf", width=25, height=15)## Warning: Removed 4 rows containing non-finite values (stat_boxplot).### For single
te_single_perc<-ggplot(subset(total_box_treat_leaf, Treatment=="20Te"), aes(x=as.factor(Leaf), y=PercentTe))+
geom_boxplot(fill="#d73027")+
theme_ines+
ylab(c("Number Te females"))+
xlab(c("Leaf Position"))+
#scale_x_discrete(labels=c())+
theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
te_single_perc
save_plot("./Plots/Te_single_dist_perc.pdf", width=7.5, height=7.5)
tu_single_perc<-ggplot(subset(total_box_treat_leaf, Treatment=="20Tu"), aes(x=as.factor(Leaf), y=PercentTu))+
geom_boxplot(fill="#fee090")+
theme_ines+
ylab(c("Number Tu females"))+
xlab(c("Leaf Position"))+
#scale_x_discrete(labels=c())+
theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
tu_single_perc
save_plot("./Plots/Tu_single_dist_perc.pdf", width=7.5, height=7.5)Te prefers leaves 3 and 4 when alone and Tu does not show a strong preference for any leaf. When Tu arrives first, there is a decrease in Te occupancy for leaves 3 and 4, with a corresponding Tu occupancy of these leaves.
Stats percentage
total_box_treat_leaf$TreatmentTe<-factor(as.factor(total_box_treat_leaf$Treatment), levels=c("20Te","10Te10Tu", "10Te10Tu48h","10Te48h10Tu","19Te1Tu" ,"19Te48h1Tu", "1Te19Tu", "1Te19Tu48h","20Tu" ))
total_box_treat_leaf$TreatmentTu<-factor(as.factor(total_box_treat_leaf$Treatment), levels=c("20Tu","10Te10Tu", "10Te10Tu48h","10Te48h10Tu","19Te1Tu" ,"19Te48h1Tu", "1Te19Tu", "1Te19Tu48h","20Te" ))
total_box_treat_leaf$Block<-as.factor(total_box_treat_leaf$Block)
total_box_treat_leaf$Leaf2<-as.factor(total_box_treat_leaf$Leaf)
str(total_box_treat_leaf)## 'data.frame': 360 obs. of 16 variables:
## $ Treatment : chr "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
## $ Box : int 1 1 1 1 2 2 2 2 3 3 ...
## $ Leaf : int 2 3 4 5 2 3 4 5 2 3 ...
## $ Block : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ sumTe : int 22 22 80 8 68 39 102 28 19 22 ...
## $ sumTu : int 27 51 27 23 30 34 54 1 24 32 ...
## $ PercentTe : num 0.1667 0.1667 0.6061 0.0606 0.2869 ...
## $ PercentTu : num 0.211 0.398 0.211 0.18 0.252 ...
## $ diffTe : num 0.0459 -0.1074 0.155 -0.0935 0.1662 ...
## $ diffTu : num 0.0436 0.1188 -0.0968 -0.0656 0.0848 ...
## $ Treatment2 : Factor w/ 7 levels "1Te19Tu","10Te10Tu",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ Treatment3 : chr "10:10\nsame time" "10:10\nsame time" "10:10\nsame time" "10:10\nsame time" ...
## $ Treatment4 : Factor w/ 7 levels "1:19\nsame time",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ TreatmentTe: Factor w/ 9 levels "20Te","10Te10Tu",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ TreatmentTu: Factor w/ 9 levels "20Tu","10Te10Tu",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ Leaf2 : Factor w/ 4 levels "2","3","4","5": 1 2 3 4 1 2 3 4 1 2 ...# adding the total mites
Total_mites<-t(sapply(c(1:length(total_box_treat_leaf[,1])), function(x){
aux<-subset(total_box_treat_leaf, Treatment==total_box_treat_leaf$Treatment[x] & Box==total_box_treat_leaf$Box[x])
te<-sum(aux$sumTe, na.rm = TRUE)
tu<-sum(aux$sumTu, na.rm = TRUE)
c(te, tu)
}))
colnames(Total_mites)<-c("TotalTe", "TotalTu")
total_box_treat_leaf<-as.data.frame(cbind(total_box_treat_leaf[, c(1:16)], Total_mites))
forTe<-subset(total_box_treat_leaf, Treatment!="20Tu")
forTe<-droplevels(forTe)
forTu<-subset(total_box_treat_leaf, Treatment!="20Te")
forTu<-droplevels(forTu)
str(forTe)## 'data.frame': 320 obs. of 18 variables:
## $ Treatment : chr "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
## $ Box : int 1 1 1 1 2 2 2 2 3 3 ...
## $ Leaf : int 2 3 4 5 2 3 4 5 2 3 ...
## $ Block : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ sumTe : int 22 22 80 8 68 39 102 28 19 22 ...
## $ sumTu : int 27 51 27 23 30 34 54 1 24 32 ...
## $ PercentTe : num 0.1667 0.1667 0.6061 0.0606 0.2869 ...
## $ PercentTu : num 0.211 0.398 0.211 0.18 0.252 ...
## $ diffTe : num 0.0459 -0.1074 0.155 -0.0935 0.1662 ...
## $ diffTu : num 0.0436 0.1188 -0.0968 -0.0656 0.0848 ...
## $ Treatment2 : Factor w/ 7 levels "1Te19Tu","10Te10Tu",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ Treatment3 : chr "10:10\nsame time" "10:10\nsame time" "10:10\nsame time" "10:10\nsame time" ...
## $ Treatment4 : Factor w/ 7 levels "1:19\nsame time",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ TreatmentTe: Factor w/ 8 levels "20Te","10Te10Tu",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ TreatmentTu: Factor w/ 8 levels "10Te10Tu","10Te10Tu48h",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Leaf2 : Factor w/ 4 levels "2","3","4","5": 1 2 3 4 1 2 3 4 1 2 ...
## $ TotalTe : int 132 132 132 132 237 237 237 237 117 117 ...
## $ TotalTu : int 128 128 128 128 119 119 119 119 133 133 ...## using cbind
# when performing the glmer the block shows variance 0, so I will run the model without the random effect
perc_te<-glm(cbind(sumTe, TotalTe)~TreatmentTe*Leaf2, data=forTe, family=binomial)
summary(perc_te)##
## Call:
## glm(formula = cbind(sumTe, TotalTe) ~ TreatmentTe * Leaf2, family = binomial,
## data = forTe)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -14.8605 -3.1028 -0.6477 1.7965 19.5147
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.126712 0.050493 -42.119 < 2e-16 ***
## TreatmentTe10Te10Tu 0.264994 0.075427 3.513 0.000443 ***
## TreatmentTe10Te10Tu48h 0.195713 0.087790 2.229 0.025793 *
## TreatmentTe10Te48h10Tu 0.175420 0.069427 2.527 0.011514 *
## TreatmentTe19Te1Tu 0.003126 0.064908 0.048 0.961594
## TreatmentTe19Te48h1Tu 0.292952 0.060735 4.823 1.41e-06 ***
## TreatmentTe1Te19Tu -0.121945 0.173747 -0.702 0.482771
## TreatmentTe1Te19Tu48h 0.845778 0.175911 4.808 1.52e-06 ***
## Leaf23 0.882468 0.061348 14.385 < 2e-16 ***
## Leaf24 1.304447 0.058644 22.243 < 2e-16 ***
## Leaf25 0.250555 0.067780 3.697 0.000219 ***
## TreatmentTe10Te10Tu:Leaf23 -0.629343 0.097135 -6.479 9.23e-11 ***
## TreatmentTe10Te10Tu48h:Leaf23 -0.096881 0.107832 -0.898 0.368951
## TreatmentTe10Te48h10Tu:Leaf23 -0.198307 0.085555 -2.318 0.020455 *
## TreatmentTe19Te1Tu:Leaf23 -0.092391 0.079242 -1.166 0.243639
## TreatmentTe19Te48h1Tu:Leaf23 -0.245901 0.074697 -3.292 0.000995 ***
## TreatmentTe1Te19Tu:Leaf23 -0.308667 0.219374 -1.407 0.159417
## TreatmentTe1Te19Tu48h:Leaf23 -1.051544 0.254258 -4.136 3.54e-05 ***
## TreatmentTe10Te10Tu:Leaf24 -0.176381 0.088773 -1.987 0.046936 *
## TreatmentTe10Te10Tu48h:Leaf24 -0.591232 0.107040 -5.523 3.32e-08 ***
## TreatmentTe10Te48h10Tu:Leaf24 -0.293750 0.081950 -3.585 0.000338 ***
## TreatmentTe19Te1Tu:Leaf24 -0.135555 0.075777 -1.789 0.073636 .
## TreatmentTe19Te48h1Tu:Leaf24 -0.453507 0.071798 -6.316 2.68e-10 ***
## TreatmentTe1Te19Tu:Leaf24 0.269364 0.197225 1.366 0.172010
## TreatmentTe1Te19Tu48h:Leaf24 -1.239908 0.242570 -5.112 3.20e-07 ***
## TreatmentTe10Te10Tu:Leaf25 -0.195058 0.103570 -1.883 0.059654 .
## TreatmentTe10Te10Tu48h:Leaf25 0.257564 0.114521 2.249 0.024508 *
## TreatmentTe10Te48h10Tu:Leaf25 0.018060 0.093096 0.194 0.846182
## TreatmentTe19Te1Tu:Leaf25 0.411658 0.084863 4.851 1.23e-06 ***
## TreatmentTe19Te48h1Tu:Leaf25 -0.225335 0.082754 -2.723 0.006470 **
## TreatmentTe1Te19Tu:Leaf25 0.378054 0.219555 1.722 0.085086 .
## TreatmentTe1Te19Tu48h:Leaf25 -0.623230 0.267245 -2.332 0.019698 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 9140.8 on 319 degrees of freedom
## Residual deviance: 5961.6 on 288 degrees of freedom
## AIC: 7590.8
##
## Number of Fisher Scoring iterations: 5perc_tu<-glm(cbind(sumTu, TotalTu)~TreatmentTu*Leaf2, data=forTu, family=binomial)
summary(perc_te)##
## Call:
## glm(formula = cbind(sumTe, TotalTe) ~ TreatmentTe * Leaf2, family = binomial,
## data = forTe)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -14.8605 -3.1028 -0.6477 1.7965 19.5147
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.126712 0.050493 -42.119 < 2e-16 ***
## TreatmentTe10Te10Tu 0.264994 0.075427 3.513 0.000443 ***
## TreatmentTe10Te10Tu48h 0.195713 0.087790 2.229 0.025793 *
## TreatmentTe10Te48h10Tu 0.175420 0.069427 2.527 0.011514 *
## TreatmentTe19Te1Tu 0.003126 0.064908 0.048 0.961594
## TreatmentTe19Te48h1Tu 0.292952 0.060735 4.823 1.41e-06 ***
## TreatmentTe1Te19Tu -0.121945 0.173747 -0.702 0.482771
## TreatmentTe1Te19Tu48h 0.845778 0.175911 4.808 1.52e-06 ***
## Leaf23 0.882468 0.061348 14.385 < 2e-16 ***
## Leaf24 1.304447 0.058644 22.243 < 2e-16 ***
## Leaf25 0.250555 0.067780 3.697 0.000219 ***
## TreatmentTe10Te10Tu:Leaf23 -0.629343 0.097135 -6.479 9.23e-11 ***
## TreatmentTe10Te10Tu48h:Leaf23 -0.096881 0.107832 -0.898 0.368951
## TreatmentTe10Te48h10Tu:Leaf23 -0.198307 0.085555 -2.318 0.020455 *
## TreatmentTe19Te1Tu:Leaf23 -0.092391 0.079242 -1.166 0.243639
## TreatmentTe19Te48h1Tu:Leaf23 -0.245901 0.074697 -3.292 0.000995 ***
## TreatmentTe1Te19Tu:Leaf23 -0.308667 0.219374 -1.407 0.159417
## TreatmentTe1Te19Tu48h:Leaf23 -1.051544 0.254258 -4.136 3.54e-05 ***
## TreatmentTe10Te10Tu:Leaf24 -0.176381 0.088773 -1.987 0.046936 *
## TreatmentTe10Te10Tu48h:Leaf24 -0.591232 0.107040 -5.523 3.32e-08 ***
## TreatmentTe10Te48h10Tu:Leaf24 -0.293750 0.081950 -3.585 0.000338 ***
## TreatmentTe19Te1Tu:Leaf24 -0.135555 0.075777 -1.789 0.073636 .
## TreatmentTe19Te48h1Tu:Leaf24 -0.453507 0.071798 -6.316 2.68e-10 ***
## TreatmentTe1Te19Tu:Leaf24 0.269364 0.197225 1.366 0.172010
## TreatmentTe1Te19Tu48h:Leaf24 -1.239908 0.242570 -5.112 3.20e-07 ***
## TreatmentTe10Te10Tu:Leaf25 -0.195058 0.103570 -1.883 0.059654 .
## TreatmentTe10Te10Tu48h:Leaf25 0.257564 0.114521 2.249 0.024508 *
## TreatmentTe10Te48h10Tu:Leaf25 0.018060 0.093096 0.194 0.846182
## TreatmentTe19Te1Tu:Leaf25 0.411658 0.084863 4.851 1.23e-06 ***
## TreatmentTe19Te48h1Tu:Leaf25 -0.225335 0.082754 -2.723 0.006470 **
## TreatmentTe1Te19Tu:Leaf25 0.378054 0.219555 1.722 0.085086 .
## TreatmentTe1Te19Tu48h:Leaf25 -0.623230 0.267245 -2.332 0.019698 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 9140.8 on 319 degrees of freedom
## Residual deviance: 5961.6 on 288 degrees of freedom
## AIC: 7590.8
##
## Number of Fisher Scoring iterations: 5Anova(perc_te)summary(perc_tu)##
## Call:
## glm(formula = cbind(sumTu, TotalTu) ~ TreatmentTu * Leaf2, family = binomial,
## data = forTu)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -16.6458 -2.2388 -0.5321 1.7716 14.2185
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.79096 0.07472 -23.970 < 2e-16 ***
## TreatmentTu10Te10Tu -0.09638 0.10333 -0.933 0.350990
## TreatmentTu10Te10Tu48h -0.32784 0.10068 -3.256 0.001129 **
## TreatmentTu10Te48h10Tu 0.05277 0.12066 0.437 0.661830
## TreatmentTu19Te1Tu -0.11858 0.27810 -0.426 0.669824
## TreatmentTu19Te48h1Tu -1.67997 0.42122 -3.988 6.65e-05 ***
## TreatmentTu1Te19Tu 0.26059 0.10530 2.475 0.013330 *
## TreatmentTu1Te19Tu48h -0.13117 0.08980 -1.461 0.144105
## Leaf23 0.57117 0.09531 5.993 2.06e-09 ***
## Leaf24 0.54654 0.09566 5.713 1.11e-08 ***
## Leaf25 0.40387 0.09786 4.127 3.67e-05 ***
## TreatmentTu10Te10Tu:Leaf23 0.06039 0.13116 0.460 0.645221
## TreatmentTu10Te10Tu48h:Leaf23 -0.02677 0.12839 -0.208 0.834854
## TreatmentTu10Te48h10Tu:Leaf23 0.16676 0.15183 1.098 0.272043
## TreatmentTu19Te1Tu:Leaf23 0.21159 0.34449 0.614 0.539074
## TreatmentTu19Te48h1Tu:Leaf23 2.22204 0.44307 5.015 5.30e-07 ***
## TreatmentTu1Te19Tu:Leaf23 -0.43990 0.13977 -3.147 0.001648 **
## TreatmentTu1Te19Tu48h:Leaf23 0.16020 0.11368 1.409 0.158770
## TreatmentTu10Te10Tu:Leaf24 -0.05358 0.13274 -0.404 0.686488
## TreatmentTu10Te10Tu48h:Leaf24 0.78653 0.12355 6.366 1.94e-10 ***
## TreatmentTu10Te48h10Tu:Leaf24 0.03358 0.15409 0.218 0.827483
## TreatmentTu19Te1Tu:Leaf24 0.39444 0.33847 1.165 0.243868
## TreatmentTu19Te48h1Tu:Leaf24 1.57372 0.45407 3.466 0.000529 ***
## TreatmentTu1Te19Tu:Leaf24 -0.14865 0.13665 -1.088 0.276694
## TreatmentTu1Te19Tu48h:Leaf24 0.26005 0.11364 2.288 0.022122 *
## TreatmentTu10Te10Tu:Leaf25 0.33045 0.13215 2.501 0.012400 *
## TreatmentTu10Te10Tu48h:Leaf25 0.18662 0.12988 1.437 0.150763
## TreatmentTu10Te48h10Tu:Leaf25 -0.61563 0.17125 -3.595 0.000325 ***
## TreatmentTu19Te1Tu:Leaf25 -0.40387 0.39127 -1.032 0.301981
## TreatmentTu19Te48h1Tu:Leaf25 1.46793 0.46070 3.186 0.001441 **
## TreatmentTu1Te19Tu:Leaf25 -0.41296 0.14362 -2.875 0.004035 **
## TreatmentTu1Te19Tu48h:Leaf25 0.01307 0.11747 0.111 0.911384
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 5399.3 on 315 degrees of freedom
## Residual deviance: 4469.2 on 284 degrees of freedom
## AIC: 5776.8
##
## Number of Fisher Scoring iterations: 5Anova(perc_tu)testInteractions(perc_te, fixed = c("Leaf2"), pairwise = c("TreatmentTe"),adjustment = "fdr")testInteractions(perc_tu, fixed = c("Leaf2"), pairwise = c("TreatmentTu"),adjustment = "fdr")Fig S4
total_box_single_24old<-subset(single_df, (Leaf==2 & leaf_age=="old") | (Leaf==4 & leaf_age=="old") | (Leaf==3 & leaf_age=="new") | (Leaf==5 & leaf_age=="new")) %>%
group_by(Treatment, Box, Leaf, Block) %>%
summarize(sumTe=sum(Te, na.rm=TRUE),sumTu=sum(Tu, na.rm=TRUE))%>%
as.data.frame()## `summarise()` has grouped output by 'Treatment', 'Box', 'Leaf'. You can override using the `.groups` argument.total_box_single_24new<-subset(single_df, (Leaf==3 & leaf_age=="old") | (Leaf==5 & leaf_age=="old") | (Leaf==2 & leaf_age=="new") | (Leaf==4 & leaf_age=="new")) %>%
group_by(Treatment, Box, Leaf, Block) %>%
summarize(sumTe=sum(Te, na.rm=TRUE),sumTu=sum(Tu, na.rm=TRUE))%>%
as.data.frame()## `summarise()` has grouped output by 'Treatment', 'Box', 'Leaf'. You can override using the `.groups` argument.# Same for the treatments
total_box_treat_24old<-subset(df_wsingle, (Leaf==2 & leaf_age=="old") | (Leaf==4 & leaf_age=="old") | (Leaf==3 & leaf_age=="new") | (Leaf==5 & leaf_age=="new")) %>%
group_by(Treatment, Box, Leaf, Block) %>%
summarize(sumTe=sum(Te, na.rm=TRUE),sumTu=sum(Tu, na.rm=TRUE))%>%
as.data.frame()## `summarise()` has grouped output by 'Treatment', 'Box', 'Leaf'. You can override using the `.groups` argument.total_box_treat_24new<-subset(df_wsingle, (Leaf==3 & leaf_age=="old") | (Leaf==5 & leaf_age=="old") | (Leaf==2 & leaf_age=="new") | (Leaf==4 & leaf_age=="new")) %>%
group_by(Treatment, Box, Leaf, Block) %>%
summarize(sumTe=sum(Te, na.rm=TRUE),sumTu=sum(Tu, na.rm=TRUE))%>%
as.data.frame()## `summarise()` has grouped output by 'Treatment', 'Box', 'Leaf'. You can override using the `.groups` argument.## Calculating percentage
#For single
total_box_single_24old$PercentTe<-sapply(c(1:length(total_box_single_24old[,1])), function(x){
aux<-subset(total_box_single_24old, Treatment==total_box_single_24old$Treatment[x] & Box==total_box_single_24old$Box[x])
aux_sum<-sum(aux$sumTe[1:4], na.rm=TRUE)
a<-total_box_single_24old$sumTe[x]/aux_sum
if(a=="NaN"){
a<-NA
}
a
})
total_box_single_24old$PercentTu<-sapply(c(1:length(total_box_single_24old[,1])), function(x){
aux<-subset(total_box_single_24old, Treatment==total_box_single_24old$Treatment[x] & Box==total_box_single_24old$Box[x])
aux_sum<-sum(aux$sumTu[1:4], na.rm=TRUE)
a<-total_box_single_24old$sumTu[x]/aux_sum
if(a=="NaN"){
a<-NA
}
a
})
total_box_single_24new$PercentTe<-sapply(c(1:length(total_box_single_24new[,1])), function(x){
aux<-subset(total_box_single_24new, Treatment==total_box_single_24new$Treatment[x] & Box==total_box_single_24new$Box[x])
aux_sum<-sum(aux$sumTe[1:4], na.rm=TRUE)
a<-total_box_single_24new$sumTe[x]/aux_sum
if(a=="NaN"){
a<-NA
}
a
})
total_box_single_24new$PercentTu<-sapply(c(1:length(total_box_single_24new[,1])), function(x){
aux<-subset(total_box_single_24new, Treatment==total_box_single_24new$Treatment[x] & Box==total_box_single_24new$Box[x])
aux_sum<-sum(aux$sumTu[1:4], na.rm=TRUE)
a<-total_box_single_24new$sumTu[x]/aux_sum
if(a=="NaN"){
a<-NA
}
a
})
# For treatments
total_box_treat_24old$PercentTe<-sapply(c(1:length(total_box_treat_24old[,1])), function(x){
aux<-subset(total_box_treat_24old, Treatment==total_box_treat_24old$Treatment[x] & Box==total_box_treat_24old$Box[x])
aux_sum<-sum(aux$sumTe[1:4], na.rm=TRUE)
a<-total_box_treat_24old$sumTe[x]/aux_sum
if(a=="NaN"){
a<-NA
}
a
})
total_box_treat_24old$PercentTu<-sapply(c(1:length(total_box_treat_24old[,1])), function(x){
aux<-subset(total_box_treat_24old, Treatment==total_box_treat_24old$Treatment[x] & Box==total_box_treat_24old$Box[x])
aux_sum<-sum(aux$sumTu[1:4], na.rm=TRUE)
a<-total_box_treat_24old$sumTu[x]/aux_sum
if(a=="NaN"){
a<-NA
}
a
})
total_box_treat_24new$PercentTe<-sapply(c(1:length(total_box_treat_24new[,1])), function(x){
aux<-subset(total_box_treat_24new, Treatment==total_box_treat_24new$Treatment[x] & Box==total_box_treat_24new$Box[x])
aux_sum<-sum(aux$sumTe[1:4], na.rm=TRUE)
a<-total_box_treat_24new$sumTe[x]/aux_sum
if(a=="NaN"){
a<-NA
}
a
})
total_box_treat_24new$PercentTu<-sapply(c(1:length(total_box_treat_24new[,1])), function(x){
aux<-subset(total_box_treat_24new, Treatment==total_box_treat_24new$Treatment[x] & Box==total_box_treat_24new$Box[x])
aux_sum<-sum(aux$sumTu[1:4], na.rm=TRUE)
a<-total_box_treat_24new$sumTu[x]/aux_sum
if(a=="NaN"){
a<-NA
}
a
})
### Now calculating the difference to the average percentage of the single regime
# first creating average data base for the single regimes
mean_percentage_single_24old<-total_box_single_24old %>%
group_by(Treatment, Leaf)%>%
summarize(meanTe=mean(PercentTe), meanTu=mean(PercentTu))%>%
as.data.frame()## `summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.mean_percentage_single_24new<-total_box_single_24new %>%
group_by(Treatment, Leaf)%>%
summarize(meanTe=mean(PercentTe), meanTu=mean(PercentTu))%>%
as.data.frame()## `summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.#Now adding to the percent data frame
aux_diffTe_24old<-t(sapply(c(1:length(total_box_treat_24old[,1])), function(x){
aux<-subset(mean_percentage_single_24old, Leaf==total_box_treat_24old$Leaf[x])
te<-total_box_treat_24old$PercentTe[x]-aux$meanTe[1]
tu<-total_box_treat_24old$PercentTu[x]-aux$meanTu[2] #because the first is always NA
c(te, tu)
}))
aux_diffTe_24new<-t(sapply(c(1:length(total_box_treat_24new[,1])), function(x){
aux<-subset(mean_percentage_single_24new, Leaf==total_box_treat_24new$Leaf[x])
te<-total_box_treat_24new$PercentTe[x]-aux$meanTe[1]
tu<-total_box_treat_24new$PercentTu[x]-aux$meanTu[2] #because the first is always NA
c(te, tu)
}))
colnames(aux_diffTe_24old)<-c("diffTe", "diffTu")
colnames(aux_diffTe_24new)<-c("diffTe", "diffTu")
total_box_treat_24old<-as.data.frame(cbind(total_box_treat_24old[,c(1:8)], aux_diffTe_24old))
str(total_box_treat_24old)## 'data.frame': 140 obs. of 10 variables:
## $ Treatment: chr "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
## $ Box : int 1 1 1 1 3 3 3 3 5 5 ...
## $ Leaf : int 2 3 4 5 2 3 4 5 2 3 ...
## $ Block : int 1 1 1 1 1 1 1 1 1 1 ...
## $ sumTe : int 22 22 80 8 19 22 44 32 77 54 ...
## $ sumTu : int 27 51 27 23 24 32 33 44 6 21 ...
## $ PercentTe: num 0.1667 0.1667 0.6061 0.0606 0.1624 ...
## $ PercentTu: num 0.211 0.398 0.211 0.18 0.18 ...
## $ diffTe : num 0.0264 -0.0446 0.109 -0.0908 0.0221 ...
## $ diffTu : num 0.14013 -0.00689 0.05903 -0.19226 0.10964 ...total_box_treat_24new<-as.data.frame(cbind(total_box_treat_24new[,c(1:8)], aux_diffTe_24new))
str(total_box_treat_24new)## 'data.frame': 140 obs. of 10 variables:
## $ Treatment: chr "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
## $ Box : int 2 2 2 2 4 4 4 4 6 6 ...
## $ Leaf : int 2 3 4 5 2 3 4 5 2 3 ...
## $ Block : int 1 1 1 1 1 1 1 1 2 2 ...
## $ sumTe : int 68 39 102 28 49 38 85 142 32 28 ...
## $ sumTu : int 30 34 54 1 34 18 45 26 20 13 ...
## $ PercentTe: num 0.287 0.165 0.43 0.118 0.156 ...
## $ PercentTu: num 0.2521 0.2857 0.4538 0.0084 0.2764 ...
## $ diffTe : num 0.1857 -0.1724 0.0254 -0.0387 0.0548 ...
## $ diffTu : num -0.0117 0.1318 -0.0098 -0.1103 0.0126 ...total_box_treat_24old$Treatment2<-factor(total_box_treat_24old$Treatment, levels=c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"))
total_box_treat_24old$Treatment3<-mapvalues(total_box_treat_24old$Treatment,c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"), c("1:19\nsame time", "10:10\nsame time", "19:1\nsame time", "1:19\nTu first", "10:10\nTu first","10:10\nTe first", "19:1\nTe first") )
total_box_treat_24old$Treatment4<-factor(total_box_treat_24old$Treatment3, levels=c("1:19\nsame time", "10:10\nsame time", "19:1\nsame time", "1:19\nTu first", "10:10\nTu first","10:10\nTe first", "19:1\nTe first") )
total_box_treat_24new$Treatment2<-factor(total_box_treat_24new$Treatment, levels=c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"))
total_box_treat_24new$Treatment3<-mapvalues(total_box_treat_24new$Treatment,c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"), c("1:19\nsame time", "10:10\nsame time", "19:1\nsame time", "1:19\nTu first", "10:10\nTu first","10:10\nTe first", "19:1\nTe first") )
total_box_treat_24new$Treatment4<-factor(total_box_treat_24new$Treatment3, levels=c("1:19\nsame time", "10:10\nsame time", "19:1\nsame time", "1:19\nTu first", "10:10\nTu first","10:10\nTe first", "19:1\nTe first") )
### For single
te_single_24old<-ggplot(subset(total_box_single_24old, Treatment=="20Te"), aes(x=as.factor(Leaf), y=PercentTe))+
geom_boxplot(fill="#d73027")+
theme_ines+
ylab(c("Number Te females"))+
xlab(c("Leaf Position"))+
#scale_x_discrete(labels=c())+
theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
te_single_24old
save_plot("./Plots/Te_single_dist_24old.pdf", width=7.5, height=7.5)
te_single_24new<-ggplot(subset(total_box_single_24new, Treatment=="20Te"), aes(x=as.factor(Leaf), y=PercentTe))+
geom_boxplot(fill="#d73027")+
theme_ines+
ylab(c("Number Te females"))+
xlab(c("Leaf Position"))+
#scale_x_discrete(labels=c())+
theme(legend.title = element_text(size=16, face="bnew"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
te_single_24new
save_plot("./Plots/Te_single_dist_24new.pdf", width=7.5, height=7.5)
tu_single_24old<-ggplot(subset(total_box_single_24old, Treatment=="20Tu"), aes(x=as.factor(Leaf), y=PercentTu))+
geom_boxplot(fill="#fee090")+
theme_ines+
ylab(c("Number Tu females"))+
xlab(c("Leaf Position"))+
ylim(c(0,1))+
#scale_x_discrete(labels=c())+
theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
tu_single_24old
save_plot("./Plots/Tu_single_dist_24old.pdf", width=7.5, height=7.5)
tu_single_24new<-ggplot(subset(total_box_single_24new, Treatment=="20Tu"), aes(x=as.factor(Leaf), y=PercentTu))+
geom_boxplot(fill="#fee090")+
theme_ines+
ylab(c("Number Tu females"))+
xlab(c("Leaf Position"))+
#scale_x_discrete(labels=c())+
ylim(c(0,1))+
theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
tu_single_24new
save_plot("./Plots/Tu_single_dist_24new.pdf", width=7.5, height=7.5)
#### Single graph with both old and new BUT WITH THE DIFFERENCES TO THE RESPECTIVE CONTROLS
new_total_box<-as.data.frame(rbind(total_box_treat_24new, total_box_treat_24new))
one_diffplot_Te<-ggplot(subset(new_total_box,Treatment!="20Te" & Treatment!="20Tu"), aes(x=as.factor(Leaf), y=diffTe, fill=Treatment4))+
facet_grid(.~Treatment4)+
geom_hline(yintercept = 0)+
geom_boxplot()+
theme_ines+
scale_fill_manual(values = c("#e0f3f8", "#91bfdb", "#4575b4", "#fee090", "#ffffbf","#d73027", "#fc8d59"), name="Treatment")+
ylab(c("Percentage of T. evansi per leaf\nrelative to single control"))+
xlab(c("Leaf Position"))+
#scale_x_discrete(labels=c())+
theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
one_diffplot_Te
save_plot("./Plots/one_diffplot_Te.pdf", width=25, height=15)
one_diffplot_Tu<-ggplot(subset(new_total_box,Treatment!="20Tu" & Treatment!="20Tu"), aes(x=as.factor(Leaf), y=diffTu, fill=Treatment4))+
facet_grid(.~Treatment4)+
geom_hline(yintercept = 0)+
geom_boxplot()+
theme_ines+
scale_fill_manual(values = c("#e0f3f8", "#91bfdb", "#4575b4", "#fee090", "#ffffbf","#d73027", "#fc8d59"), name="Treatment")+
ylab(c("Percentage of T. urticae per leaf\nrelative to single control"))+
xlab(c("Leaf Position"))+
#scale_x_discreTu(labels=c())+
theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
one_diffplot_Tu
save_plot("./Plots/one_diffplot_Tu.pdf", width=25, height=15)Statistics (accounting for leaf order)
Distribution
hist(new_total_box$diffTe)
hist(new_total_box$diffTu)
descdist(new_total_box$diffTe, discrete = FALSE, boot=500)
## summary statistics
## ------
## min: -0.4049489 max: 0.5098335
## median: -0.01485167
## mean: 6.147063e-18
## estimated sd: 0.1639964
## estimated skewness: 0.4349104
## estimated kurtosis: 3.35588descdist(new_total_box$diffTu, discrete = FALSE, boot=500)
## summary statistics
## ------
## min: -0.463581 max: 0.8461161
## median: -0.03990551
## mean: -6.938894e-18
## estimated sd: 0.2168845
## estimated skewness: 1.133723
## estimated kurtosis: 5.350042#Log normal fits bothStatistics
str(new_total_box)## 'data.frame': 280 obs. of 13 variables:
## $ Treatment : chr "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
## $ Box : int 2 2 2 2 4 4 4 4 6 6 ...
## $ Leaf : int 2 3 4 5 2 3 4 5 2 3 ...
## $ Block : int 1 1 1 1 1 1 1 1 2 2 ...
## $ sumTe : int 68 39 102 28 49 38 85 142 32 28 ...
## $ sumTu : int 30 34 54 1 34 18 45 26 20 13 ...
## $ PercentTe : num 0.287 0.165 0.43 0.118 0.156 ...
## $ PercentTu : num 0.2521 0.2857 0.4538 0.0084 0.2764 ...
## $ diffTe : num 0.1857 -0.1724 0.0254 -0.0387 0.0548 ...
## $ diffTu : num -0.0117 0.1318 -0.0098 -0.1103 0.0126 ...
## $ Treatment2: Factor w/ 7 levels "1Te19Tu","10Te10Tu",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ Treatment3: chr "10:10\nsame time" "10:10\nsame time" "10:10\nsame time" "10:10\nsame time" ...
## $ Treatment4: Factor w/ 7 levels "1:19\nsame time",..: 2 2 2 2 2 2 2 2 2 2 ...# Making leaf and Block as factors
new_total_box$Leaf2<-as.factor(new_total_box$Leaf)
new_total_box$Block2<-as.factor(new_total_box$Block)
mla_Te<-lmer(log(diffTe+1)~Treatment2*Leaf2+ (1|Block2), data=new_total_box)## boundary (singular) fit: see ?isSingularmla_Tu<-lmer(log(diffTu+1)~Treatment2*Leaf2+ (1|Block2), data=new_total_box)## boundary (singular) fit: see ?isSingularsummary(mla_Te)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(diffTe + 1) ~ Treatment2 * Leaf2 + (1 | Block2)
## Data: new_total_box
##
## REML criterion at convergence: -172.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7410 -0.5897 -0.1344 0.5338 3.0181
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block2 (Intercept) 0.00000 0.0000
## Residual 0.02289 0.1513
## Number of obs: 280, groups: Block2, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.04177 0.04784 252.00000 -0.873 0.38343
## Treatment210Te10Tu 0.11006 0.06766 252.00000 1.627 0.10506
## Treatment219Te1Tu 0.08407 0.06766 252.00000 1.243 0.21518
## Treatment21Te19Tu48h 0.12896 0.06766 252.00000 1.906 0.05780 .
## Treatment210Te10Tu48h 0.12091 0.06766 252.00000 1.787 0.07514 .
## Treatment210Te48h10Tu 0.13183 0.06766 252.00000 1.948 0.05248 .
## Treatment219Te48h1Tu 0.16474 0.06766 252.00000 2.435 0.01560 *
## Leaf23 -0.06911 0.06766 252.00000 -1.021 0.30800
## Leaf24 0.08866 0.06766 252.00000 1.310 0.19124
## Leaf25 0.05220 0.06766 252.00000 0.772 0.44111
## Treatment210Te10Tu:Leaf23 -0.17565 0.09569 252.00000 -1.836 0.06759 .
## Treatment219Te1Tu:Leaf23 -0.00421 0.09569 252.00000 -0.044 0.96494
## Treatment21Te19Tu48h:Leaf23 -0.13432 0.09569 252.00000 -1.404 0.16164
## Treatment210Te10Tu48h:Leaf23 0.03848 0.09569 252.00000 0.402 0.68791
## Treatment210Te48h10Tu:Leaf23 -0.06526 0.09569 252.00000 -0.682 0.49586
## Treatment219Te48h1Tu:Leaf23 -0.08879 0.09569 252.00000 -0.928 0.35436
## Treatment210Te10Tu:Leaf24 -0.10078 0.09569 252.00000 -1.053 0.29326
## Treatment219Te1Tu:Leaf24 -0.05353 0.09569 252.00000 -0.559 0.57635
## Treatment21Te19Tu48h:Leaf24 -0.28975 0.09569 252.00000 -3.028 0.00272 **
## Treatment210Te10Tu48h:Leaf24 -0.25875 0.09569 252.00000 -2.704 0.00731 **
## Treatment210Te48h10Tu:Leaf24 -0.27732 0.09569 252.00000 -2.898 0.00408 **
## Treatment219Te48h1Tu:Leaf24 -0.27320 0.09569 252.00000 -2.855 0.00466 **
## Treatment210Te10Tu:Leaf25 -0.11585 0.09569 252.00000 -1.211 0.22716
## Treatment219Te1Tu:Leaf25 -0.21379 0.09569 252.00000 -2.234 0.02634 *
## Treatment21Te19Tu48h:Leaf25 -0.09856 0.09569 252.00000 -1.030 0.30396
## Treatment210Te10Tu48h:Leaf25 -0.20139 0.09569 252.00000 -2.105 0.03631 *
## Treatment210Te48h10Tu:Leaf25 -0.12585 0.09569 252.00000 -1.315 0.18962
## Treatment219Te48h1Tu:Leaf25 -0.22977 0.09569 252.00000 -2.401 0.01706 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 28 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularsummary(mla_Tu)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(diffTu + 1) ~ Treatment2 * Leaf2 + (1 | Block2)
## Data: new_total_box
##
## REML criterion at convergence: -59.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3319 -0.3926 -0.0802 0.3748 3.1119
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block2 (Intercept) 0.00000 0.0000
## Residual 0.03577 0.1891
## Number of obs: 280, groups: Block2, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.812e-02 5.981e-02 2.520e+02 0.637 0.524468
## Treatment210Te10Tu -1.184e-01 8.458e-02 2.520e+02 -1.400 0.162797
## Treatment219Te1Tu -7.500e-02 8.458e-02 2.520e+02 -0.887 0.376065
## Treatment21Te19Tu48h -1.472e-01 8.458e-02 2.520e+02 -1.741 0.082929
## Treatment210Te10Tu48h -1.154e-01 8.458e-02 2.520e+02 -1.364 0.173674
## Treatment210Te48h10Tu 4.201e-02 8.458e-02 2.520e+02 0.497 0.619878
## Treatment219Te48h1Tu -3.390e-01 8.458e-02 2.520e+02 -4.008 8.06e-05
## Leaf23 4.902e-04 8.458e-02 2.520e+02 0.006 0.995380
## Leaf24 -2.032e-01 8.458e-02 2.520e+02 -2.402 0.017026
## Leaf25 -4.833e-03 8.458e-02 2.520e+02 -0.057 0.954481
## Treatment210Te10Tu:Leaf23 1.201e-01 1.196e-01 2.520e+02 1.004 0.316407
## Treatment219Te1Tu:Leaf23 -5.645e-02 1.196e-01 2.520e+02 -0.472 0.637413
## Treatment21Te19Tu48h:Leaf23 1.064e-01 1.196e-01 2.520e+02 0.890 0.374532
## Treatment210Te10Tu48h:Leaf23 6.027e-02 1.196e-01 2.520e+02 0.504 0.614793
## Treatment210Te48h10Tu:Leaf23 -3.288e-02 1.196e-01 2.520e+02 -0.275 0.783627
## Treatment219Te48h1Tu:Leaf23 4.796e-01 1.196e-01 2.520e+02 4.010 8.02e-05
## Treatment210Te10Tu:Leaf24 1.025e-01 1.196e-01 2.520e+02 0.857 0.392153
## Treatment219Te1Tu:Leaf24 2.287e-01 1.196e-01 2.520e+02 1.912 0.057062
## Treatment21Te19Tu48h:Leaf24 3.930e-01 1.196e-01 2.520e+02 3.286 0.001162
## Treatment210Te10Tu48h:Leaf24 4.161e-01 1.196e-01 2.520e+02 3.479 0.000594
## Treatment210Te48h10Tu:Leaf24 1.753e-02 1.196e-01 2.520e+02 0.147 0.883606
## Treatment219Te48h1Tu:Leaf24 5.114e-01 1.196e-01 2.520e+02 4.276 2.71e-05
## Treatment210Te10Tu:Leaf25 2.398e-01 1.196e-01 2.520e+02 2.005 0.046033
## Treatment219Te1Tu:Leaf25 1.919e-02 1.196e-01 2.520e+02 0.160 0.872649
## Treatment21Te19Tu48h:Leaf25 1.129e-01 1.196e-01 2.520e+02 0.944 0.346168
## Treatment210Te10Tu48h:Leaf25 -7.036e-03 1.196e-01 2.520e+02 -0.059 0.953143
## Treatment210Te48h10Tu:Leaf25 -1.267e-01 1.196e-01 2.520e+02 -1.059 0.290719
## Treatment219Te48h1Tu:Leaf25 2.014e-01 1.196e-01 2.520e+02 1.684 0.093459
##
## (Intercept)
## Treatment210Te10Tu
## Treatment219Te1Tu
## Treatment21Te19Tu48h .
## Treatment210Te10Tu48h
## Treatment210Te48h10Tu
## Treatment219Te48h1Tu ***
## Leaf23
## Leaf24 *
## Leaf25
## Treatment210Te10Tu:Leaf23
## Treatment219Te1Tu:Leaf23
## Treatment21Te19Tu48h:Leaf23
## Treatment210Te10Tu48h:Leaf23
## Treatment210Te48h10Tu:Leaf23
## Treatment219Te48h1Tu:Leaf23 ***
## Treatment210Te10Tu:Leaf24
## Treatment219Te1Tu:Leaf24 .
## Treatment21Te19Tu48h:Leaf24 **
## Treatment210Te10Tu48h:Leaf24 ***
## Treatment210Te48h10Tu:Leaf24
## Treatment219Te48h1Tu:Leaf24 ***
## Treatment210Te10Tu:Leaf25 *
## Treatment219Te1Tu:Leaf25
## Treatment21Te19Tu48h:Leaf25
## Treatment210Te10Tu48h:Leaf25
## Treatment210Te48h10Tu:Leaf25
## Treatment219Te48h1Tu:Leaf25 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 28 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularAnova(mla_Te)Anova(mla_Tu)#### Contrasts
levels(as.factor(new_total_box$Treatment))## [1] "10Te10Tu" "10Te10Tu48h" "10Te48h10Tu" "19Te1Tu" "19Te48h1Tu"
## [6] "1Te19Tu" "1Te19Tu48h"mat_Leaf2<-matrix(ncol=10, nrow=7, 0)
#Compare Tu arriving first 10:10
mat_Leaf2[1,1]<-1
mat_Leaf2[2,1]<- -1
#Compare Te arriving first 10:10
mat_Leaf2[1,2]<-1
mat_Leaf2[3,2]<- -1
#Compare Te arriving first 19:1
mat_Leaf2[4,3]<-1
mat_Leaf2[5,3]<- -1
#Compare Tu arriving first 1:19
mat_Leaf2[6,4]<-1
mat_Leaf2[7,4]<- -1
#Compare Same arrival density effect Te
mat_Leaf2[1,5]<-1
mat_Leaf2[4,5]<- -1
#Compare Same arrival density effect Tu
mat_Leaf2[1,6]<-1
mat_Leaf2[6,6]<- -1
#Compare Tu arriving first (both densities)
mat_Leaf2[1,7]<-1
mat_Leaf2[2,7]<- -1
mat_Leaf2[6,7]<-1
mat_Leaf2[7,7]<- -1
#Compare Te arriving first (both densities)
mat_Leaf2[4,8]<-1
mat_Leaf2[5,8]<- -1
mat_Leaf2[1,8]<-1
mat_Leaf2[3,8]<- -1
#Compare 1:19Tu (both timepoints)
mat_Leaf2[1,9]<-1
mat_Leaf2[2,9]<- 1
mat_Leaf2[6,9]<- -1
mat_Leaf2[7,9]<- -1
#Compare 19Te:1 (both timepoints)
mat_Leaf2[1,10]<-1
mat_Leaf2[3,10]<- 1
mat_Leaf2[4,10]<- -1
mat_Leaf2[5,10]<- -1
rownames(mat_Leaf2)<-c(levels(as.factor(df_wsingle$Treatment)))
colnames(mat_Leaf2)<-c("Tufirst_10:10", "Tefirst_10:10", "Tefirst_19:1", "Tufirst_1:19", "DensityTe19", "DensityTu19", "Tu arriving first", "Te arriving first", "Tu density", "Te density")
mat_Leaf2## Tufirst_10:10 Tefirst_10:10 Tefirst_19:1 Tufirst_1:19 DensityTe19
## 10Te10Tu 1 1 0 0 1
## 10Te10Tu48h -1 0 0 0 0
## 10Te48h10Tu 0 -1 0 0 0
## 19Te1Tu 0 0 1 0 -1
## 19Te48h1Tu 0 0 -1 0 0
## 1Te19Tu 0 0 0 1 0
## 1Te19Tu48h 0 0 0 -1 0
## DensityTu19 Tu arriving first Te arriving first Tu density
## 10Te10Tu 1 1 1 1
## 10Te10Tu48h 0 -1 0 1
## 10Te48h10Tu 0 0 -1 0
## 19Te1Tu 0 0 1 0
## 19Te48h1Tu 0 0 -1 0
## 1Te19Tu -1 1 0 -1
## 1Te19Tu48h 0 -1 0 -1
## Te density
## 10Te10Tu 1
## 10Te10Tu48h 0
## 10Te48h10Tu 1
## 19Te1Tu -1
## 19Te48h1Tu -1
## 1Te19Tu 0
## 1Te19Tu48h 0test_la_Te<-testInteractions(mla_Te, fixed=c("Leaf2"), custom = list(Treatment2=mat_Leaf2), adjustment = "fdr" )
test_la_Tu<-testInteractions(mla_Tu, fixed=c("Leaf2"), custom = list(Treatment2=mat_Leaf2), adjustment = "fdr" )
test_la_Tetest_la_TuCalculating agreggation scores to test if changes between leaf, leaf age and direction
We will use C-score (package bipartite, which normalises c_score between 0 and 1) to test for aggregation between the two species for the different treatments
C score
For the C score we need a presence absence matrix per leaflet inside of a leaf
str(df)## 'data.frame': 3214 obs. of 17 variables:
## $ Box : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : chr "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
## $ Block : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Leaf : int 2 2 2 2 2 2 3 3 3 3 ...
## $ Leaflet : int 1 2 3 1 2 3 1 2 3 4 ...
## $ Direction : Factor w/ 2 levels "Down","Up": 1 1 1 2 2 2 1 1 1 1 ...
## $ Te : int 4 2 1 6 6 3 0 0 4 0 ...
## $ Tu : int 1 2 5 1 12 6 1 2 0 2 ...
## $ leaf_age : chr "old" "old" "old" "old" ...
## $ Density : chr "10:10" "10:10" "10:10" "10:10" ...
## $ Timing : chr "same time" "same time" "same time" "same time" ...
## $ ratio : num 4 1 0.2 6 0.5 0.5 1 0.5 4 0.5 ...
## $ Total : int 5 4 6 7 18 9 1 2 4 2 ...
## $ proportionTe: num 0.8 0.5 0.167 0.857 0.333 ...
## $ proportionTu: num 0.2 0.5 0.833 0.143 0.667 ...
## $ Block2 : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ Box2 : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...df$Leaf_leaflet<-sapply(c(1:length(df[,1])), function(x){paste(df$Leaf[x],df$Leaflet[x], sep="_")})
df$Leaf_leaflet_dir<-sapply(c(1:length(df[,1])), function(x){paste(df$Leaf[x],df$Leaflet[x], df$Direction[x], sep="_")})
### Summing everything per leaflet
aux_df2<-df %>%
group_by(Treatment, Density, Timing , Block, Box, Leaf, Leaflet) %>%
summarise(Te=sum(Te, na.rm=TRUE), Tu=sum(Tu, na.rm=TRUE)) %>%
as.data.frame()## `summarise()` has grouped output by 'Treatment', 'Density', 'Timing', 'Block', 'Box', 'Leaf'. You can override using the `.groups` argument.aux_df2$Leaf_leaflet<-sapply(c(1:length(aux_df2[,1])), function(x){paste(aux_df2$Leaf[x],aux_df2$Leaflet[x], sep="_")})
aux_df2$Tu[which(aux_df2$Tu >1)]<-1
aux_df2$Te[which(aux_df2$Te >1)]<-1
str(aux_df2)## 'data.frame': 1607 obs. of 10 variables:
## $ Treatment : chr "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
## $ Density : chr "10:10" "10:10" "10:10" "10:10" ...
## $ Timing : chr "same time" "same time" "same time" "same time" ...
## $ Block : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Box : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Leaf : int 2 2 2 3 3 3 3 3 4 4 ...
## $ Leaflet : int 1 2 3 1 2 3 4 5 1 2 ...
## $ Te : num 1 1 1 0 0 1 1 1 1 1 ...
## $ Tu : num 1 1 1 1 1 1 1 1 1 1 ...
## $ Leaf_leaflet: chr "2_1" "2_2" "2_3" "3_1" ...tu_df2<-pivot_wider(aux_df2[,c(1, 5, 9,10)], names_from=Leaf_leaflet, values_from=Tu)
te_df2<-pivot_wider(aux_df2[,c(1, 5, 8,10)], names_from=Leaf_leaflet, values_from=Te)
tu_df2$species<-"Tu"
te_df2$species<-"Te"
mat_02<-as.data.frame(rbind(tu_df2, te_df2))
mat_02mat_02[is.na(mat_02)]<-0
c_score_mat<-expand.grid(Box=c(1:10), Treatment=levels(as.factor(df$Treatment)))
c_score_mat$cscore<-NA
str(mat_02)## 'data.frame': 180 obs. of 21 variables:
## $ Treatment: chr "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
## $ Box : int 1 2 3 4 5 6 7 8 9 10 ...
## $ 2_1 : num 1 1 1 1 1 1 1 1 0 1 ...
## $ 2_2 : num 1 1 1 1 1 1 1 1 1 1 ...
## $ 2_3 : num 1 1 1 1 1 1 1 1 1 1 ...
## $ 3_1 : num 1 0 1 0 0 0 1 1 1 1 ...
## $ 3_2 : num 1 1 0 1 1 0 1 1 1 1 ...
## $ 3_3 : num 1 1 1 1 1 1 1 1 1 1 ...
## $ 3_4 : num 1 1 1 1 1 1 1 0 1 1 ...
## $ 3_5 : num 1 1 0 1 1 1 1 0 1 1 ...
## $ 4_1 : num 1 1 1 1 1 1 1 1 1 0 ...
## $ 4_2 : num 1 1 1 1 0 1 1 1 1 1 ...
## $ 4_3 : num 1 1 1 1 1 1 1 1 0 1 ...
## $ 4_4 : num 1 1 0 1 0 1 1 1 0 0 ...
## $ 4_5 : num 1 1 1 1 1 1 1 1 1 1 ...
## $ 5_1 : num 1 0 1 1 1 1 1 1 1 1 ...
## $ 5_2 : num 1 0 1 1 0 1 1 1 1 1 ...
## $ 5_3 : num 1 0 1 1 1 1 1 1 1 1 ...
## $ 5_4 : num 1 0 1 1 0 1 1 1 1 1 ...
## $ 5_5 : num 1 1 1 1 1 1 1 1 1 1 ...
## $ species : chr "Tu" "Tu" "Tu" "Tu" ...c_score_mat$cscore<-sapply(c(1:length(c_score_mat$Box)), function(x){
C.score(t(subset(mat_02, Box==c_score_mat$Box[x] & as.character(Treatment)==as.character(c_score_mat$Treatment[x]))[,c(3:20)]), normalise = TRUE)
})
c_score_mat$Treatment2<-factor(c_score_mat$Treatment, levels=c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"))
c_score_mat$Treatment3<-mapvalues(c_score_mat$Treatment2,c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"), c("1:19\nsame time", "10:10\nsame time", "19:1\nsame time", "1:19\nTu first", "10:10\nTu first","10:10\nTe first", "19:1\nTe first") )
c_score_mat<-c_score_mat[-which(is.na(c_score_mat)),]
c_score_mat<-c_score_mat[-which(c_score_mat$cscore=="NaN"),]
c_score_mat<-droplevels(c_score_mat)Figure S5
ggplot(c_score_mat, aes(x=Treatment3, y=cscore, fill=Treatment3))+
geom_boxplot(na.rm=TRUE)+
theme_ines+
scale_fill_manual(values = c("#e0f3f8", "#91bfdb", "#4575b4", "#fee090", "#ffffbf","#d73027", "#fc8d59"), name="Treatment")+
scale_x_discrete(labels=c("","","","","","",""))+
ylab("C-Score")+
xlab("Treatment")+
theme(axis.text.x = element_text(angle=60, vjust = 0.5))
save_plot("./Plots/C_score.pdf")Test if C scores differ between treatments
shapiro.test(c_score_mat$cscore)##
## Shapiro-Wilk normality test
##
## data: c_score_mat$cscore
## W = 0.62001, p-value = 4.254e-11hist(c_score_mat$cscore)
descdist(c_score_mat$cscore+1, discrete = FALSE, boot=500)
## summary statistics
## ------
## min: 1 max: 2
## median: 1
## mean: 1.166668
## estimated sd: 0.29382
## estimated skewness: 2.014817
## estimated kurtosis: 6.210139plot(fitdist(c_score_mat$cscore+1, "gamma"))## $start.arg
## $start.arg$shape
## [1] 16.03821
##
## $start.arg$rate
## [1] 13.74702
##
##
## $fix.arg
## NULL
plot(fitdist(c_score_mat$cscore, "beta", "mme"))## $start.arg
## $start.arg$shape1
## [1] 0.1060946
##
## $start.arg$shape2
## [1] 0.5304664
##
##
## $fix.arg
## NULL
plot(fitdist(c_score_mat$cscore+1, "lnorm"))## $start.arg
## $start.arg$meanlog
## [1] 0.1294786
##
## $start.arg$sdlog
## [1] 0.2097881
##
##
## $fix.arg
## NULL
fit_dist_cscore<-summary(fitdist(c_score_mat$cscore, "beta", "mme"))## $start.arg
## $start.arg$shape1
## [1] 0.1060946
##
## $start.arg$shape2
## [1] 0.5304664
##
##
## $fix.arg
## NULLstr(fit_dist_cscore)## List of 20
## $ estimate : Named num [1:2] 0.106 0.53
## ..- attr(*, "names")= chr [1:2] "shape1" "shape2"
## $ method : chr "mme"
## $ sd : NULL
## $ cor : NULL
## $ vcov : NULL
## $ loglik : num Inf
## $ aic : num -Inf
## $ bic : num -Inf
## $ n : int 59
## $ data : num [1:59] 1 0 0 0.333 0.444 ...
## $ distname : chr "beta"
## $ fix.arg : NULL
## $ fix.arg.fun: NULL
## $ dots : NULL
## $ convergence: num 0
## $ discrete : logi FALSE
## $ weights : NULL
## $ ddistname : chr "dbeta"
## $ pdistname : chr "pbeta"
## $ qdistname : chr "qbeta"
## - attr(*, "class")= chr [1:2] "summary.fitdist" "fitdist"c_score_mat$Block<-mapvalues(c_score_mat$Box, c(1,2,3,4,5,6,7,8,9,10), c(1,1,1,1,1,2,2,2,2,2))
c_score_mat$Block<-as.factor(c_score_mat$Block)
cs2<-glm((cscore+1)~Treatment, family=Gamma, data=c_score_mat)
summary(cs2)##
## Call:
## glm(formula = (cscore + 1) ~ Treatment, family = Gamma, data = c_score_mat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.35779 -0.13385 -0.03956 0.00000 0.56236
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.68354 0.06092 11.220 1.68e-15 ***
## Treatment10Te10Tu48h 0.14278 0.08560 1.668 0.10135
## Treatment10Te48h10Tu 0.10912 0.08390 1.301 0.19915
## Treatment19Te1Tu 0.31646 0.09834 3.218 0.00222 **
## Treatment19Te48h1Tu 0.25675 0.09865 2.603 0.01203 *
## Treatment1Te19Tu 0.15264 0.08393 1.819 0.07473 .
## Treatment1Te19Tu48h 0.27742 0.09007 3.080 0.00331 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.04766633)
##
## Null deviance: 2.9115 on 58 degrees of freedom
## Residual deviance: 2.1399 on 52 degrees of freedom
## AIC: 3.3808
##
## Number of Fisher Scoring iterations: 4Anova(cs2, type=2)testInteractions(cs2, pairwise="Treatment", adjustment = "fdr")Order of arrival affects probability of coexistence?
Calculating FD from Godoy and Levine 2014
The demographic parameter in our case is lambdai -1 / lambdaj -1 and the competitive response the alphas which are arranged in a different way to calculate niche diff.
aij describes the per capita effect of species j on species i –> so in the script, in the results the row is the focal so the Te inter is the effect that Tu has on Te, and Tu inter vice versa.
fitness ratio = ((nj-1)/(ni-1))* sqrt(aij/ajj * aii/aji) Niche overlap = sqrt(aij/ajj* aji/aii)
Don’t forget that niche differences is 1- niche overlap
The data is generated by the script alpha_estimate_fixed_final.Rmd, so please run this script before hand
df_alpha_lambda<-read.csv("./Estimate_Lambda_fixed_allTreats/Results_all.csv", header=TRUE)
ggplot(df_alpha_lambda, aes(x=Comp, y=Alpha, colour=Arrival, fill=Arrival))+
facet_grid(.~Species)+
geom_errorbar(aes(ymin=Alpha_Lower, ymax=Alpha_Upper), position=position_dodge(0.5), size=0.75)+
geom_point(size=2, position=position_dodge(0.5), shape=21, colour="black")+
geom_line(aes(group=Arrival), position=position_dodge2(0.5))+
scale_colour_manual(values = c("#4575b4","#d73027", "#fee090"))+
scale_fill_manual(values =c("#4575b4","#d73027", "#fee090"))+
theme_ines+
ylab("Competitive ability")+
xlab("")+
theme(legend.position = "bottom")
save_plot("./Estimate_Lambda_fixed_allTreats/alpha.pdf")
ggplot(df_alpha_lambda, aes(x=Species, y=Lambda, colour=Arrival, fill=Arrival))+
geom_hline(yintercept =2, colour="lightgray")+
geom_hline(yintercept =3.5, colour="lightgray")+
geom_errorbar(aes(ymin=Lambda_Lower, ymax=Lambda_Upper), position=position_dodge(0.5), size=0.75)+
geom_point(size=2, position=position_dodge(0.5), shape=21, colour="black")+
scale_colour_manual(values = c("#4575b4","#d73027", "#fee090"))+
scale_fill_manual(values = c("#4575b4","#d73027", "#fee090"))+
theme_ines+
ylab("Lambda Growth Rate")+
xlab("")+
theme(legend.position = "bottom")
save_plot("./Estimate_Lambda_fixed_allTreats/lambda.pdf")
coex_df2<-as.data.frame(matrix(nrow=3, ncol=7))
colnames(coex_df2)<-c("Treatment", "NO", "NO_L", "NO_U", "FD", "FD_L","FD_U")
coex_df2$Treatment<-c("same time","Te first","Tu first")
str(coex_df2)## 'data.frame': 3 obs. of 7 variables:
## $ Treatment: chr "same time" "Te first" "Tu first"
## $ NO : logi NA NA NA
## $ NO_L : logi NA NA NA
## $ NO_U : logi NA NA NA
## $ FD : logi NA NA NA
## $ FD_L : logi NA NA NA
## $ FD_U : logi NA NA NAstr(df_alpha_lambda)## 'data.frame': 12 obs. of 9 variables:
## $ Species : chr "Te" "Tu" "Te" "Tu" ...
## $ Arrival : chr "same time" "same time" "Te first" "Te first" ...
## $ Comp : chr "Inter" "Inter" "Inter" "Inter" ...
## $ Lambda : num 3.6 1.73 3.6 1.73 3.6 ...
## $ Lambda_Upper: num 4.18 1.9 4.18 1.9 4.18 ...
## $ Lambda_Lower: num 3.02 1.56 3.02 1.56 3.02 ...
## $ Alpha : num 0.0532 0.0261 0.0586 0.0455 0.0679 ...
## $ Alpha_Upper : num 0.0648 0.0368 0.0693 0.0731 0.0767 ...
## $ Alpha_Lower : num 0.0417 0.0153 0.048 0.0178 0.0591 ...head(df_alpha_lambda)### Important!
#The order here is Te inter (aji), Tu inter (aij), Te intra (ajj), Tu intra (aii) for each treatment
coex_df2$NO<-sapply(c(1:3), function(x){
aux<-subset(df_alpha_lambda, as.character(Arrival)==as.character(coex_df2$Treatment[x]))
#Assuming that the order is always te tu, te tu, inter inter, intra intra
nd<-(sqrt(aux$Alpha[2]/aux$Alpha[3]* aux$Alpha[1]/aux$Alpha[4]))
nd
})
coex_df2$NO_L<-sapply(c(1:3), function(x){
aux<-subset(df_alpha_lambda, as.character(Arrival)==as.character(coex_df2$Treatment[x]))
#Assuming that the order is always te tu, te tu, inter inter, intra intra
nd<- sqrt(aux$Alpha_Lower[2]/aux$Alpha_Lower[3]* aux$Alpha_Lower[1]/aux$Alpha_Lower[4])
nd
})
coex_df2$NO_U<-sapply(c(1:3), function(x){
aux<-subset(df_alpha_lambda, as.character(Arrival)==as.character(coex_df2$Treatment[x]))
#Assuming that the order is always te tu, te tu, inter inter, intra intra
nd<-(sqrt(aux$Alpha_Upper[2]/aux$Alpha_Upper[3]* aux$Alpha_Upper[1]/aux$Alpha_Upper[4]))
nd
})
#fd= ((nj-1)/(ni-1))* sqrt(aij/ajj * aii/aji)
coex_df2$FD<-sapply(c(1:3), function(x){
aux<-subset(df_alpha_lambda, as.character(Arrival)==as.character(coex_df2$Treatment[x]))
#The order here is Te inter (aji), Tu inter (aij), Te intra (ajj), Tu intra (aii)
fd<-((aux$Lambda[1]-1)/(aux$Lambda[2]-1))*(sqrt(aux$Alpha[2]/aux$Alpha[3])* sqrt(aux$Alpha[4]/aux$Alpha[1]))
fd
})
coex_df2$FD_L<-sapply(c(1:3), function(x){
aux<-subset(df_alpha_lambda, as.character(Arrival)==as.character(coex_df2$Treatment[x]))
#Assuming that the order is always te tu, te tu, inter inter, intra intra
fd<-((aux$Lambda_Lower[1]-1)/(aux$Lambda_Lower[2]-1))*(sqrt(aux$Alpha_Lower[2]/aux$Alpha_Lower[3])* sqrt(aux$Alpha_Lower[4]/aux$Alpha_Lower[1]))
fd
})
coex_df2$FD_U<-sapply(c(1:3), function(x){
aux<-subset(df_alpha_lambda, as.character(Arrival)==as.character(coex_df2$Treatment[x]))
#Assuming that the order is always te tu, te tu, inter inter, intra intra
fd<-((aux$Lambda_Upper[1]-1)/(aux$Lambda_Upper[2]-1))*(sqrt(aux$Alpha_Upper[2]/aux$Alpha_Upper[3])* sqrt(aux$Alpha_Upper[4]/aux$Alpha_Upper[1]))
fd
})
coex_df2$ND<-1-coex_df2$NO
coex_df2$ND_L<-1-coex_df2$NO_L
coex_df2$ND_U<-1-coex_df2$NO_U
bound_df<-data.frame(niche_overlap=c(seq(0,3, 0.05))) # creating a vector with niche overlap
bound_df$niche_diff<-(1-bound_df$niche_overlap) # calculating stabilizating differences from niche overlap 1-rho
bound_df$fitness_differences_sp_1<-(1/bound_df$niche_overlap) # solid line in your graph this is ok
bound_df$fitness_differences_sp_temp<- 1-bound_df$fitness_differences_sp_1 #this is an intermediate step to see the differences above one
#which is later incorporated into the 2 species
bound_df$fitness_differences_sp_2<- 1+ bound_df$fitness_differences_sp_temp# dashed line in your graph which is NOT OK in your original graph
#Here I added the difference calculated in the previous step to make the curves simetric. they must be simmetric!
#otherwise one species has more chance of priority effects than the other.
bound_df<-bound_df[, -4]
#remove the intermediate step
str(bound_df)## 'data.frame': 61 obs. of 4 variables:
## $ niche_overlap : num 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 ...
## $ niche_diff : num 1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 ...
## $ fitness_differences_sp_1: num Inf 20 10 6.67 5 ...
## $ fitness_differences_sp_2: num -Inf -18 -8 -4.67 -3 ...colnames(bound_df)<-c("NO", "ND", "FD", "FD2")Fig 2 - Coexistence plot
ggplot(coex_df2,aes(x=ND, y=FD, colour=Treatment) )+
geom_hline(yintercept = 1 ,colour="lightgray")+
geom_vline(xintercept = 0 , colour="lightgray")+
geom_ribbon(data=bound_df, aes(x=ND, ymin=FD,ymax=FD2), colour="Black", fill="lightgrey", alpha=0.5)+
geom_line(data=bound_df, aes(x=ND, y=FD), colour="black")+
geom_line(data=bound_df, aes(x=ND, y=FD2),colour="black")+
geom_errorbar(data=coex_df2, aes(ymin=FD_L, ymax=FD_U), colour="black", width=0.1)+
geom_errorbarh(data=coex_df2, aes(xmin=ND_L, xmax=ND_U), colour="black", height=0.1)+
geom_point(data=coex_df2, size=2.75, shape=21, aes(fill=Treatment), colour="black")+
ylab(c(""))+
xlab(c(""))+
#scale_x_discrete(labels=c())+
scale_colour_manual(values = c("#4575b4","#d73027", "#fee090"), name="")+
scale_fill_manual(values = c("#4575b4","#d73027", "#fee090"), name="")+
scale_x_continuous(breaks = seq(-1,1, 0.2), expand=expansion(c(0,0)), limits = c(-1.5, 1))+
ylim(c(-1.5,5))+
theme_bw()+
theme_ines+
theme(legend.position = c(0.9,0.9), legend.background = element_rect(fill=NA), legend.text = element_text(size=10), axis.text = element_text(colour="black"))+
geom_richtext(label="<span style='color:black'>Priority effects", x=-0.82, y=1.1, size=4, color=NA, fill=NA)+
geom_richtext(label="<span style='color:black'>*T. urticae* excluded", x=0, y=2.5, size=4, color=NA, fill=NA)+
geom_richtext(label="<span style='color:black'>*T. evansi* excluded", x=0, y=0, size=4, color=NA, fill=NA)+
geom_richtext(label="<span style='color:black'>Coexistence", x=0.56, y=1.1, size=4, color=NA, fill=NA)+
expand_limits(x = c(-1, 0.75), y = 0)+
coord_cartesian(ylim = c(0, 5), xlim = c(-1,0.7))## Warning: Removed 10 row(s) containing missing values (geom_path).
## Warning: Removed 10 row(s) containing missing values (geom_path).
save_plot("./Estimate_Lambda_fixed_allTreats/coex.pdf")## Warning: Removed 10 row(s) containing missing values (geom_path).
## Warning: Removed 10 row(s) containing missing values (geom_path).